GlobalTrackingTools/src/DynamicTruncation.cc

0001 /**
0002  *  Class: DynamicTruncation
0003  *
0004  *  Description:
0005  *  class for the dynamical stop of the KF according to the
0006  *  compatibility degree between the extrapolated track
0007  *  state and the reconstructed segment in the muon chambers
0008  *
0009  *  Authors :
0010  *  D. Pagano & G. Bruno - UCL Louvain
0011  *
0012  *  \modified by C. Caputo, UCLouvain
0013  **/
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 #include "RecoMuon/GlobalTrackingTools/interface/DynamicTruncation.h"
0017 #include "DataFormats/MuonDetId/interface/RPCDetId.h"
0018 #include "TrackingTools/TrackFitters/interface/RecHitLessByDet.h"
0019 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0020 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
0021 #include "RecoMuon/Records/interface/MuonRecoGeometryRecord.h"
0022 #include "RecoMuon/Navigation/interface/MuonNavigationSchool.h"
0023 #include "RecoMuon/Navigation/interface/MuonNavigationPrinter.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0025 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
0026 #include "RecoMuon/TrackingTools/interface/MuonPatternRecoDumper.h"
0027 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0028 #include "TrackPropagation/SteppingHelixPropagator/interface/SteppingHelixPropagator.h"
0029
0030 #define MAX_THR 1e7
0031
0032 using namespace edm;
0033 using namespace std;
0034 using namespace reco;
0035
0036 namespace dyt_utils {
0037   static const std::map<etaRegion, std::string> etaRegionStr{{etaRegion::eta0p8, "eta0p8"},
0038                                                              {etaRegion::eta1p2, "eta1p2"},
0039                                                              {etaRegion::eta2p0, "eta2p0"},
0040                                                              {etaRegion::eta2p2, "eta2p2"},
0041                                                              {etaRegion::eta2p4, "eta2p4"}};
0042 };
0043
0044 DynamicTruncation::Config::Config(edm::ConsumesCollector iC)
0045     : cscGeomToken_(iC.esConsumes()),
0046       muonRecHitBuilderToken_(iC.esConsumes(edm::ESInputTag("", "MuonRecHitBuilder"))),
0047       updatorToken_(iC.esConsumes(edm::ESInputTag("", "KFUpdator"))),
0048       navMuonToken_(iC.esConsumes()),
0049       dytThresholdsToken_(iC.esConsumes()),
0050       dtAlignmentErrorsToken_(iC.esConsumes()),
0051       cscAlignmentErrorsToken_(iC.esConsumes()) {}
0052
0053 DynamicTruncation::DynamicTruncation(Config const &config,
0054                                      const edm::EventSetup &eventSetup,
0055                                      const MuonServiceProxy &theService) {
0056   propagator = theService.propagator("SmartPropagatorAny");
0057   propagatorPF = theService.propagator("SmartPropagatorAny");
0058   propagatorCompatibleDet = theService.propagator("SmartPropagatorAny");
0059   theG = theService.trackingGeometry();
0060   theMuonRecHitBuilder = eventSetup.getHandle(config.muonRecHitBuilderToken_);
0061   updatorHandle = eventSetup.getHandle(config.updatorToken_);
0062   cscGeom = eventSetup.getHandle(config.cscGeomToken_);
0063   navMuon = eventSetup.getHandle(config.navMuonToken_);
0064   magfield = theService.magneticField();
0065   navigation = std::make_unique<DirectMuonNavigation>(theService.detLayerGeometry());
0066   getSegs = std::make_unique<ChamberSegmentUtility>();
0067   {
0068     edm::ESHandle<DYTThrObject> dytThresholdsH = eventSetup.getHandle(config.dytThresholdsToken_);
0069     AlignmentErrorsExtended const &dtAlignmentErrorsExtended = eventSetup.getData(config.dtAlignmentErrorsToken_);
0070     AlignmentErrorsExtended const &cscAlignmentErrorsExtended = eventSetup.getData(config.cscAlignmentErrorsToken_);
0071
0072     thrManager = std::make_unique<ThrParameters>(dytThresholdsH, dtAlignmentErrorsExtended, cscAlignmentErrorsExtended);
0073   }
0074   useDBforThr = thrManager->isValidThdDB();
0075   dytThresholds = nullptr;
0076   if (useDBforThr)
0077     dytThresholds = thrManager->getInitialThresholds();
0078
0079   doUpdateOfKFStates = true;
0080   useParametrizedThr = false;
0081 }
0082
0083 DynamicTruncation::~DynamicTruncation() {}
0084
0085 void DynamicTruncation::update(TrajectoryStateOnSurface &tsos, ConstRecHitPointer rechit) {
0086   TrajectoryStateOnSurface temp = updatorHandle->update(tsos, *rechit);
0087   if (temp.isValid())
0088     tsos = updatorHandle->update(tsos, *rechit);
0089 }
0090
0091 void DynamicTruncation::updateWithDThits(TrajectoryStateOnSurface &tsos, DTRecSegment4D const &bestDTSeg) {
0092   ConstRecHitContainer tmprecHits;
0093   vector<const TrackingRecHit *> DTrh = bestDTSeg.recHits();
0094   for (vector<const TrackingRecHit *>::iterator it = DTrh.begin(); it != DTrh.end(); it++) {
0095     tmprecHits.push_back(theMuonRecHitBuilder->build(*it));
0096   }
0097   sort(tmprecHits);
0098   for (ConstRecHitContainer::const_iterator it = tmprecHits.begin(); it != tmprecHits.end(); ++it) {
0099     DTLayerId layid((*it)->det()->geographicalId());
0100     TrajectoryStateOnSurface temp = propagator->propagate(tsos, theG->idToDet(layid)->surface());
0101     if (temp.isValid()) {
0102       TrajectoryStateOnSurface tempTsos = updatorHandle->update(temp, **it);
0103       if (tempTsos.isValid())
0104         tsos = tempTsos;
0105     }
0106   }
0107 }
0108
0109 void DynamicTruncation::updateWithCSChits(TrajectoryStateOnSurface &tsos, CSCSegment const &bestCSCSeg) {
0110   ConstRecHitContainer tmprecHits;
0111   vector<CSCRecHit2D> CSCrh = bestCSCSeg.specificRecHits();
0112   for (vector<CSCRecHit2D>::iterator it = CSCrh.begin(); it != CSCrh.end(); ++it) {
0113     tmprecHits.push_back(theMuonRecHitBuilder->build(&*it));
0114   }
0115   sort(tmprecHits);
0116   for (ConstRecHitContainer::const_iterator it = tmprecHits.begin(); it != tmprecHits.end(); ++it) {
0117     const CSCLayer *cscLayer = cscGeom->layer((*it)->det()->geographicalId());
0118     TrajectoryStateOnSurface temp = propagator->propagate(tsos, cscLayer->surface());
0119     if (temp.isValid()) {
0120       TrajectoryStateOnSurface tempTsos = updatorHandle->update(temp, **it);
0121       if (tempTsos.isValid())
0122         tsos = tempTsos;
0123     }
0124   }
0125 }
0126
0127 /////////////////////////////////
0128 ///// Configuration methods /////
0129 /////////////////////////////////
0130 void DynamicTruncation::setSelector(int selector) {
0131   if (selector < 0 || selector > 2)
0132     throw cms::Exception("NotAvailable") << "DYT selector: wrong option!" << endl;
0133   //if (selector == 0) cout << "[DYT disabled]\n";
0134   //if (selector == 1) cout << "[use all compatible stations]\n";
0135   //if (selector == 2) cout << "[stop at second consecutive incompatible station]\n";
0136   DYTselector = selector;
0137 }
0138
0139 void DynamicTruncation::setUseAPE(bool useAPE_) { useAPE = useAPE_; }
0140
0141 void DynamicTruncation::setUpdateState(bool upState) { doUpdateOfKFStates = upState; }
0142
0143 void DynamicTruncation::setThr(const vector<int> &thr) {
0144   if (thr.size() == 2) {
0145     for (unsigned int i = 0; i < thr.size(); i++)
0146       if (thr[i] >= 0)
0147         Thrs.push_back(thr[i]);
0148       else
0149         Thrs.push_back(MAX_THR);
0150     return;
0151   }
0152   throw cms::Exception("NotAvailable")
0153       << "WARNING: wrong size for the threshold vector!\nExpected size: 2\n   Found size: " << thr.size();
0154 }
0155
0156 void DynamicTruncation::setThrsMap(const edm::ParameterSet &par) {
0157   for (auto const &region : dyt_utils::etaRegionStr) {
0158     parameters[region.first] = par.getParameter<std::vector<double> >(region.second);
0159   }
0160 }
0161 /////////////////////////////////
0162 /////////////////////////////////
0163 /////////////////////////////////
0164
0165 //===> filter
0166 TransientTrackingRecHit::ConstRecHitContainer DynamicTruncation::filter(const Trajectory &traj) {
0167   result.clear();
0168   prelFitMeas.clear();
0169
0170   // Get APE maps
0171   dtApeMap = thrManager->GetDTApeMap();
0172   cscApeMap = thrManager->GetCSCApeMap();
0173
0174   // Get Last tracker TSOS (updated)
0175   vector<TrajectoryMeasurement> muonMeasurements = traj.measurements();
0176   TrajectoryMeasurement lastTKm = muonMeasurements.front();
0177   for (vector<TrajectoryMeasurement>::const_iterator imT = muonMeasurements.begin(); imT != muonMeasurements.end();
0178        imT++) {
0179     if (!(*imT).recHit()->isValid())
0180       continue;
0181     const TransientTrackingRecHit *hit = &(*(*imT).recHit());
0182     if (hit->geographicalId().det() == DetId::Tracker) {
0183       result.push_back((*imT).recHit());
0184       if (!(*imT).forwardPredictedState().isValid())
0185         continue;
0186       if ((*imT).forwardPredictedState().globalPosition().mag() >
0187           lastTKm.forwardPredictedState().globalPosition().mag())
0188         lastTKm = *imT;
0189     }
0190   }
0191   currentState = lastTKm.forwardPredictedState();
0192   update(currentState, lastTKm.recHit());
0193
0194   prelFitState = lastTKm.forwardPredictedState();
0195   update(prelFitState, lastTKm.recHit());
0196   prelFitMeas = result;
0197
0198   // Run the DYT
0199   filteringAlgo();
0200
0201   return result;
0202 }
0203
0204 //===> filteringAlgo
0205 void DynamicTruncation::filteringAlgo() {
0206   map<int, vector<DetId> > compatibleIds;
0207   map<int, vector<DTRecSegment4D> > dtSegMap;
0208   map<int, vector<CSCSegment> > cscSegMap;
0209   int incompConLay = 0;
0210   nStationsUsed = 0;
0211
0212   // Get list of compatible layers
0213   compatibleDets(currentState, compatibleIds);
0214
0215   // Fill segment maps
0216   fillSegmentMaps(compatibleIds, dtSegMap, cscSegMap);
0217
0218   // Do a preliminary fit
0219   if (useDBforThr)
0220     preliminaryFit(compatibleIds, dtSegMap, cscSegMap);
0221
0222   // Loop on compatible layers
0223   for (map<int, vector<DetId> >::iterator it = compatibleIds.begin(); it != compatibleIds.end(); ++it) {
0224     int stLayer = stationfromDet(it->second.front());
0225     DTRecSegment4D bestDTSeg;
0226     CSCSegment bestCSCSeg;
0227     double bestDTEstimator = MAX_THR;
0228     double bestCSCEstimator = MAX_THR;
0229     vector<DTRecSegment4D> dtSegs = dtSegMap[it->first];
0230     vector<CSCSegment> cscSegs = cscSegMap[it->first];
0231
0232     // DT case: find the most compatible segment
0233     TrajectoryStateOnSurface tsosDTlayer;
0234     testDTstation(currentState, dtSegs, bestDTEstimator, bestDTSeg, tsosDTlayer);
0235
0236     // CSC case: find the most compatible segment
0237     TrajectoryStateOnSurface tsosCSClayer;
0238     testCSCstation(currentState, cscSegs, bestCSCEstimator, bestCSCSeg, tsosCSClayer);
0239
0240     // Decide whether to keep the layer or not
0241     bool chosenLayer =
0242         chooseLayers(incompConLay, bestDTEstimator, bestDTSeg, tsosDTlayer, bestCSCEstimator, bestCSCSeg, tsosCSClayer);
0243     fillDYTInfos(stLayer, chosenLayer, incompConLay, bestDTEstimator, bestCSCEstimator, bestDTSeg, bestCSCSeg);
0244   }
0245   //cout << "Number of used stations = " << nStationsUsed << endl;
0246 }
0247
0248 //===> stationfromDet
0249 int DynamicTruncation::stationfromDet(DetId const &det) {
0250   if (det.subdetId() == MuonSubdetId::CSC) {
0251     CSCDetId ch(det);
0252     return ch.station();
0253   }
0254   if (det.subdetId() == MuonSubdetId::DT) {
0255     DTChamberId ch(det);
0256     return ch.station();
0257   }
0258   return 0;
0259 }
0260
0261 //===> fillDYTInfos
0262 void DynamicTruncation::fillDYTInfos(int const &st,
0263                                      bool const &chosenLayer,
0264                                      int &incompConLay,
0265                                      double const &bestDTEstimator,
0266                                      double const &bestCSCEstimator,
0267                                      DTRecSegment4D const &bestDTSeg,
0268                                      CSCSegment const &bestCSCSeg) {
0269   if (chosenLayer) {
0270     nStationsUsed++;
0271     incompConLay = 0;
0272     if (bestDTEstimator <= bestCSCEstimator) {
0273       estimatorMap[st] = bestDTEstimator;
0274       DetId id(bestDTSeg.chamberId());
0275       idChamberMap[st] = id;
0276     } else {
0277       DetId id(bestCSCSeg.cscDetId());
0278       idChamberMap[st] = id;
0279       estimatorMap[st] = bestCSCEstimator;
0280     }
0281     usedStationMap[st] = true;
0282   } else {
0283     incompConLay++;
0284     estimatorMap[st] = -1;
0285     usedStationMap[st] = false;
0286   }
0287 }
0288
0289 //===> compatibleDets
0290 void DynamicTruncation::compatibleDets(TrajectoryStateOnSurface &tsos, map<int, vector<DetId> > &detMap) {
0291   MuonPatternRecoDumper dumper;
0292   MeasurementEstimator *theEstimator = new Chi2MeasurementEstimator(1000, 1000);
0293   vector<const DetLayer *> navLayers;
0294   navLayers = navigation->compatibleLayers(*(currentState.freeState()), alongMomentum);
0295   unsigned int ilayerCorrected = 0;
0296   for (unsigned int ilayer = 0; ilayer < navLayers.size(); ilayer++) {
0297     // Skip RPC layers
0298     if (navLayers[ilayer]->subDetector() != GeomDetEnumerators::DT &&
0299         navLayers[ilayer]->subDetector() != GeomDetEnumerators::CSC)
0300       continue;
0301     ilayerCorrected++;
0302     vector<DetLayer::DetWithState> comps =
0303         navLayers[ilayer]->compatibleDets(currentState, *propagatorCompatibleDet, *theEstimator);
0304     //cout << comps.size() << " compatible Dets with " << navLayers[ilayer]->subDetector() << " Layer " << ilayer << " "
0305     //<< dumper.dumpLayer(navLayers[ilayer]);
0306     if (!comps.empty()) {
0307       for (unsigned int icomp = 0; icomp < comps.size(); icomp++) {
0308         DetId id(comps[icomp].first->geographicalId().rawId());
0309         detMap[ilayerCorrected].push_back(id);
0310       }
0311     }
0312   }
0313   if (theEstimator)
0314     delete theEstimator;
0315 }
0316
0317 //===> fillSegmentMaps
0318 void DynamicTruncation::fillSegmentMaps(map<int, vector<DetId> > &compatibleIds,
0319                                         map<int, vector<DTRecSegment4D> > &dtSegMap,
0320                                         map<int, vector<CSCSegment> > &cscSegMap) {
0321   for (map<int, vector<DetId> >::iterator it = compatibleIds.begin(); it != compatibleIds.end(); ++it) {
0322     vector<DetId> ids = compatibleIds[it->first];
0323     for (unsigned j = 0; j < ids.size(); j++) {
0324       if (ids[j].subdetId() == MuonSubdetId::CSC) {
0325         CSCDetId ch(ids[j]);
0326         vector<CSCSegment> tmp = getSegs->getCSCSegmentsInChamber(ch);
0327         for (unsigned int k = 0; k < tmp.size(); k++)
0328           cscSegMap[it->first].push_back(tmp[k]);
0329       }
0330       if (ids[j].subdetId() == MuonSubdetId::DT) {
0331         DTChamberId ch(ids[j]);
0332         vector<DTRecSegment4D> tmp = getSegs->getDTSegmentsInChamber(ch);
0333         for (unsigned int k = 0; k < tmp.size(); k++)
0334           dtSegMap[it->first].push_back(tmp[k]);
0335       }
0336     }
0337   }
0338 }
0339
0340 //===> testDTstation
0341 void DynamicTruncation::testDTstation(TrajectoryStateOnSurface &startingState,
0342                                       vector<DTRecSegment4D> const &segments,
0343                                       double &bestEstimator,
0344                                       DTRecSegment4D &bestSeg,
0345                                       TrajectoryStateOnSurface &tsosdt) {
0346   if (segments.empty())
0347     return;
0348   for (unsigned int iSeg = 0; iSeg < segments.size(); iSeg++) {
0349     DTChamberId chamber(segments[iSeg].chamberId());
0350     if (!propagator->propagate(startingState, theG->idToDet(chamber)->surface()).isValid())
0351       continue;
0352     tsosdt = propagator->propagate(startingState, theG->idToDet(chamber)->surface());
0353     //if (!tsosdt.isValid()) continue;
0354     LocalError apeLoc;
0355     if (useAPE)
0356       apeLoc = ErrorFrameTransformer().transform(dtApeMap.find(chamber)->second, theG->idToDet(chamber)->surface());
0357     StateSegmentMatcher estim(tsosdt, segments[iSeg], apeLoc);
0358     double estimator = estim.value();
0359     //cout << "estimator DT = " << estimator << endl;
0360     if (estimator >= bestEstimator)
0361       continue;
0362     bestEstimator = estimator;
0363     bestSeg = segments[iSeg];
0364   }
0365 }
0366
0367 //===> testCSCstation
0368 void DynamicTruncation::testCSCstation(TrajectoryStateOnSurface &startingState,
0369                                        vector<CSCSegment> const &segments,
0370                                        double &bestEstimator,
0371                                        CSCSegment &bestSeg,
0372                                        TrajectoryStateOnSurface &tsoscsc) {
0373   if (segments.empty())
0374     return;
0375   for (unsigned int iSeg = 0; iSeg < segments.size(); iSeg++) {
0376     CSCDetId chamber(segments[iSeg].cscDetId());
0377     if (!propagator->propagate(startingState, theG->idToDet(chamber)->surface()).isValid())
0378       continue;
0379     tsoscsc = propagator->propagate(startingState, theG->idToDet(chamber)->surface());
0380     //if (!tsoscsc.isValid()) continue;
0381     LocalError apeLoc;
0382     if (useAPE)
0383       apeLoc = ErrorFrameTransformer().transform(cscApeMap.find(chamber)->second, theG->idToDet(chamber)->surface());
0384     StateSegmentMatcher estim(tsoscsc, segments[iSeg], apeLoc);
0385     double estimator = estim.value();
0386     //cout << "estimator CSC = " << estimator << endl;
0387     if (estimator >= bestEstimator)
0388       continue;
0389     bestEstimator = estimator;
0390     bestSeg = segments[iSeg];
0391   }
0392 }
0393
0394 //===> useSegment
0395 void DynamicTruncation::useSegment(DTRecSegment4D const &bestDTSeg, TrajectoryStateOnSurface const &tsosDT) {
0396   result.push_back(theMuonRecHitBuilder->build(&bestDTSeg));
0397   if (doUpdateOfKFStates)
0398     updateWithDThits(currentState, bestDTSeg);
0399   else
0400     currentState = tsosDT;
0401 }
0402
0403 //===> useSegment
0404 void DynamicTruncation::useSegment(CSCSegment const &bestCSCSeg, TrajectoryStateOnSurface const &tsosCSC) {
0405   result.push_back(theMuonRecHitBuilder->build(&bestCSCSeg));
0406   if (doUpdateOfKFStates)
0407     updateWithCSChits(currentState, bestCSCSeg);
0408   else
0409     currentState = tsosCSC;
0410 }
0411
0412 //===> preliminaryFit
0413 void DynamicTruncation::preliminaryFit(map<int, vector<DetId> > compatibleIds,
0414                                        map<int, vector<DTRecSegment4D> > dtSegMap,
0415                                        map<int, vector<CSCSegment> > cscSegMap) {
0416   for (map<int, vector<DetId> >::iterator it = compatibleIds.begin(); it != compatibleIds.end(); ++it) {
0417     DTRecSegment4D bestDTSeg;
0418     CSCSegment bestCSCSeg;
0419     double bestDTEstimator = MAX_THR;
0420     double bestCSCEstimator = MAX_THR;
0421     double initThr = MAX_THR;
0422     vector<DTRecSegment4D> dtSegs = dtSegMap[it->first];
0423     vector<CSCSegment> cscSegs = cscSegMap[it->first];
0424
0425     // DT case: find the most compatible segment
0426     TrajectoryStateOnSurface tsosDTlayer;
0427     testDTstation(prelFitState, dtSegs, bestDTEstimator, bestDTSeg, tsosDTlayer);
0428
0429     // CSC case: find the most compatible segment
0430     TrajectoryStateOnSurface tsosCSClayer;
0431     testCSCstation(prelFitState, cscSegs, bestCSCEstimator, bestCSCSeg, tsosCSClayer);
0432
0433     // Decide whether to keep the layer or not
0434     if (bestDTEstimator == MAX_THR && bestCSCEstimator == MAX_THR)
0435       continue;
0436     if (bestDTEstimator <= bestCSCEstimator) {
0437       getThresholdFromCFG(initThr, DetId(bestDTSeg.chamberId()));
0438       if (bestDTEstimator >= initThr)
0439         continue;
0440       prelFitMeas.push_back(theMuonRecHitBuilder->build(&bestDTSeg));
0441       auto aSegRH = prelFitMeas.back();
0442       auto uRes = updatorHandle->update(tsosDTlayer, *aSegRH);
0443       if (uRes.isValid()) {
0444         prelFitState = uRes;
0445       } else {
0446         prelFitMeas.pop_back();
0447       }
0448     } else {
0449       getThresholdFromCFG(initThr, DetId(bestCSCSeg.cscDetId()));
0450       if (bestCSCEstimator >= initThr)
0451         continue;
0452       prelFitMeas.push_back(theMuonRecHitBuilder->build(&bestCSCSeg));
0453       auto aSegRH = prelFitMeas.back();
0454       auto uRes = updatorHandle->update(tsosCSClayer, *aSegRH);
0455       if (uRes.isValid()) {
0456         prelFitState = uRes;
0457       } else {
0458         prelFitMeas.pop_back();
0459       }
0460     }
0461   }
0462   if (!prelFitMeas.empty())
0463     prelFitMeas.pop_back();
0464   for (auto imrh = prelFitMeas.rbegin(); imrh != prelFitMeas.rend(); ++imrh) {
0465     DetId id = (*imrh)->geographicalId();
0466     TrajectoryStateOnSurface tmp = propagatorPF->propagate(prelFitState, theG->idToDet(id)->surface());
0467     if (tmp.isValid())
0468       prelFitState = tmp;
0469   }
0470   muonPTest = prelFitState.globalMomentum().perp();
0471   muonETAest = prelFitState.globalMomentum().eta();
0472 }
0473
0474 //===> chooseLayers
0475 bool DynamicTruncation::chooseLayers(int &incompLayers,
0476                                      double const &bestDTEstimator,
0477                                      DTRecSegment4D const &bestDTSeg,
0478                                      TrajectoryStateOnSurface const &tsosDT,
0479                                      double const &bestCSCEstimator,
0480                                      CSCSegment const &bestCSCSeg,
0481                                      TrajectoryStateOnSurface const &tsosCSC) {
0482   double initThr = MAX_THR;
0483   if (bestDTEstimator == MAX_THR && bestCSCEstimator == MAX_THR)
0484     return false;
0485   if (bestDTEstimator <= bestCSCEstimator) {
0486     // Get threshold for the chamber
0487     if (useDBforThr)
0488       getThresholdFromDB(initThr, DetId(bestDTSeg.chamberId()));
0489     else
0490       getThresholdFromCFG(initThr, DetId(bestDTSeg.chamberId()));
0491     if (DYTselector == 0 || (DYTselector == 1 && bestDTEstimator < initThr) ||
0492         (DYTselector == 2 && incompLayers < 2 && bestDTEstimator < initThr)) {
0493       useSegment(bestDTSeg, tsosDT);
0494       return true;
0495     }
0496   } else {
0497     // Get threshold for the chamber
0498     if (useDBforThr)
0499       getThresholdFromDB(initThr, DetId(bestCSCSeg.cscDetId()));
0500     else
0501       getThresholdFromCFG(initThr, DetId(bestCSCSeg.cscDetId()));
0502     if (DYTselector == 0 || (DYTselector == 1 && bestCSCEstimator < initThr) ||
0503         (DYTselector == 2 && incompLayers < 2 && bestCSCEstimator < initThr)) {
0504       useSegment(bestCSCSeg, tsosCSC);
0505       return true;
0506     }
0507   }
0508   return false;
0509 }
0510
0511 //===> getThresholdFromDB
0512 void DynamicTruncation::getThresholdFromDB(double &thr, DetId const &id) {
0513   vector<DYTThrObject::DytThrStruct> thrvector = dytThresholds->thrsVec;
0514   for (vector<DYTThrObject::DytThrStruct>::const_iterator it = thrvector.begin(); it != thrvector.end(); it++) {
0515     DYTThrObject::DytThrStruct obj = (*it);
0516     if (obj.id == id) {
0517       thr = obj.thr;
0518       break;
0519     }
0520   }
0521   if (useParametrizedThr)
0522     correctThrByPAndEta(thr);
0523 }
0524
0525 //===> correctThrByPAndEta
0526 void DynamicTruncation::correctThrByPAndEta(double &thr) {
0527   auto parametricThreshold = [this] {
0528     double thr50 = this->parameters[this->region].at(0);
0529     double p0 = this->parameters[this->region].at(1);
0530     double p1 = this->parameters[this->region].at(2);
0531     return thr50 * (1 + p0 * p_reco + std::pow(this->p_reco, p1));
0532   };
0533
0534   std::set<dyt_utils::etaRegion> regionsToExclude = {
0535       dyt_utils::etaRegion::eta2p0, dyt_utils::etaRegion::eta2p2, dyt_utils::etaRegion::eta2p4};
0536
0537   if (!regionsToExclude.count(this->region))
0538     thr = parametricThreshold();
0539 }
0540
0541 void DynamicTruncation::setEtaRegion() {
0542   float absEta = std::abs(eta_reco);
0543   // Defaul value for muons with abs(eta) > 2.4
0544   region = dyt_utils::etaRegion::eta2p4;
0545
0546   if (absEta <= 0.8)
0547     region = dyt_utils::etaRegion::eta0p8;
0548   else if (absEta <= 1.2)
0549     region = dyt_utils::etaRegion::eta1p2;
0550   else if (absEta <= 2.0)
0551     region = dyt_utils::etaRegion::eta2p0;
0552   else if (absEta <= 2.2)
0553     region = dyt_utils::etaRegion::eta2p2;
0554   else if (absEta <= 2.4)
0555     region = dyt_utils::etaRegion::eta2p4;
0556 }
0557
0558 //===> getThresholdFromCFG
0559 void DynamicTruncation::getThresholdFromCFG(double &thr, DetId const &id) {
0560   if (id.subdetId() == MuonSubdetId::DT) {
0561     thr = Thrs[0];
0562   }
0563   if (id.subdetId() == MuonSubdetId::CSC) {
0564     thr = Thrs[1];
0565   }
0566
0567   if (useParametrizedThr)
0568     correctThrByPAndEta(thr);
0569 }
0570
0571 //===> sort
0572 void DynamicTruncation::sort(ConstRecHitContainer &recHits) {
0573   unsigned int i = 0;
0574   unsigned int j = 0;
0575   ConstRecHitContainer::size_type n = recHits.size();
0576   for (i = 1; i < n; ++i)
0577     for (j = n - 1; j >= i; --j)
0578       if (recHits[j - 1]->globalPosition().mag() > recHits[j]->globalPosition().mag())
0579         swap(recHits[j - 1], recHits[j]);
0580 }