Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:13 

0001 /**

0002  * \file ME0SegAlgRU.cc

0003  *

0004  *  \author M. Maggi for ME0

0005  *  \from V.Palichik & N.Voytishin

0006  *  \some functions and structure taken from SK algo by M.Sani and SegFit class by T.Cox

0007  */
0008 #include "ME0SegAlgoRU.h"
0009 #include "MuonSegFit.h"
0010 #include "Geometry/GEMGeometry/interface/ME0Layer.h"
0011 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0012 #include "DataFormats/Math/interface/deltaPhi.h"
0013 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 
0016 #include <algorithm>
0017 #include <cmath>
0018 #include <iostream>
0019 #include <string>
0020 
0021 #include <Math/Functions.h>
0022 #include <Math/SVector.h>
0023 #include <Math/SMatrix.h>
0024 
0025 ME0SegAlgoRU::ME0SegAlgoRU(const edm::ParameterSet& ps) : ME0SegmentAlgorithmBase(ps), myName("ME0SegAlgoRU") {
0026   allowWideSegments = ps.getParameter<bool>("allowWideSegments");
0027   doCollisions = ps.getParameter<bool>("doCollisions");
0028 
0029   stdParameters.maxChi2Additional = ps.getParameter<double>("maxChi2Additional");
0030   stdParameters.maxChi2Prune = ps.getParameter<double>("maxChi2Prune");
0031   stdParameters.maxChi2GoodSeg = ps.getParameter<double>("maxChi2GoodSeg");
0032   stdParameters.maxPhiSeeds = ps.getParameter<double>("maxPhiSeeds");
0033   stdParameters.maxPhiAdditional = ps.getParameter<double>("maxPhiAdditional");
0034   stdParameters.maxETASeeds = ps.getParameter<double>("maxETASeeds");
0035   stdParameters.maxTOFDiff = ps.getParameter<double>("maxTOFDiff");
0036   stdParameters.requireCentralBX = ps.getParameter<bool>("requireCentralBX");
0037   stdParameters.minNumberOfHits = ps.getParameter<unsigned int>("minNumberOfHits");
0038   stdParameters.requireBeamConstr = true;
0039 
0040   wideParameters = stdParameters;
0041   wideParameters.maxChi2Prune *= 2;
0042   wideParameters.maxChi2GoodSeg *= 2;
0043   wideParameters.maxPhiSeeds *= 2;
0044   wideParameters.maxPhiAdditional *= 2;
0045   wideParameters.minNumberOfHits += 1;
0046 
0047   displacedParameters = stdParameters;
0048   displacedParameters.maxChi2Additional *= 2;
0049   displacedParameters.maxChi2Prune = 100;
0050   displacedParameters.maxChi2GoodSeg *= 5;
0051   displacedParameters.maxPhiSeeds *= 2;
0052   displacedParameters.maxPhiAdditional *= 2;
0053   displacedParameters.maxETASeeds *= 2;
0054   displacedParameters.requireBeamConstr = false;
0055 
0056   LogDebug("ME0SegAlgoRU") << myName << " has algorithm cuts set to: \n"
0057                            << "--------------------------------------------------------------------\n"
0058                            << "allowWideSegments   = " << allowWideSegments << "\n"
0059                            << "doCollisions        = " << doCollisions << "\n"
0060                            << "maxChi2Additional   = " << stdParameters.maxChi2Additional << "\n"
0061                            << "maxChi2Prune        = " << stdParameters.maxChi2Prune << "\n"
0062                            << "maxChi2GoodSeg      = " << stdParameters.maxChi2GoodSeg << "\n"
0063                            << "maxPhiSeeds         = " << stdParameters.maxPhiSeeds << "\n"
0064                            << "maxPhiAdditional    = " << stdParameters.maxPhiAdditional << "\n"
0065                            << "maxETASeeds         = " << stdParameters.maxETASeeds << "\n"
0066                            << "maxTOFDiff          = " << stdParameters.maxTOFDiff << "\n"
0067                            << std::endl;
0068 
0069   theChamber = nullptr;
0070 }
0071 
0072 std::vector<ME0Segment> ME0SegAlgoRU::run(const ME0Chamber* chamber, const HitAndPositionContainer& rechits) {
0073 #ifdef EDM_ML_DEBUG  // have lines below only compiled when in debug mode

0074   ME0DetId chId(chamber->id());
0075   edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegmentAlgorithm::run] build segments in chamber " << chId
0076                                    << " which contains " << rechits.size() << " rechits";
0077   for (unsigned int iH = 0; iH < rechits.size(); ++iH) {
0078     auto me0id = rechits[iH].rh->me0Id();
0079     auto rhLP = rechits[iH].lp;
0080     edm::LogVerbatim("ME0SegAlgoRU") << "[RecHit :: Loc x = " << std::showpos << std::setw(9) << rhLP.x()
0081                                      << " Glb y = " << std::showpos << std::setw(9) << rhLP.y()
0082                                      << " Time = " << std::showpos << rechits[iH].rh->tof() << " -- " << me0id.rawId()
0083                                      << " = " << me0id << " ]" << std::endl;
0084   }
0085 #endif
0086 
0087   if (rechits.size() < 3 || rechits.size() > 300) {
0088     return std::vector<ME0Segment>();
0089   }
0090 
0091   theChamber = chamber;
0092 
0093   std::vector<unsigned int> recHits_per_layer(6, 0);
0094   for (unsigned int iH = 0; iH < rechits.size(); ++iH) {
0095     recHits_per_layer[rechits[iH].layer - 1]++;
0096   }
0097 
0098   BoolContainer used(rechits.size(), false);
0099 
0100   // We have at least 2 hits. We intend to try all possible pairs of hits to start

0101   // segment building. 'All possible' means each hit lies on different layers in the chamber.

0102   // after all same size segs are build we get rid of the overcrossed segments using the chi2 criteria

0103   // the hits from the segs that are left are marked as used and are not added to segs in future iterations

0104   // the hits from 3p segs are marked as used separately in order to try to assamble them in longer segments

0105   // in case there is a second pass

0106 
0107   // Choose first hit (as close to IP as possible) h1 and second hit

0108   // (as far from IP as possible) h2 To do this we iterate over hits

0109   // in the chamber by layer - pick two layers.  Then we

0110   // iterate over hits within each of these layers and pick h1 and h2

0111   // these.  If they are 'close enough' we build an empty

0112   // segment.  Then try adding hits to this segment.

0113 
0114   std::vector<ME0Segment> segments;
0115 
0116   auto doStd = [&]() {
0117     for (unsigned int n_seg_min = 6u; n_seg_min >= stdParameters.minNumberOfHits; --n_seg_min)
0118       lookForSegments(stdParameters, n_seg_min, rechits, recHits_per_layer, used, segments);
0119   };
0120   auto doDisplaced = [&]() {
0121     for (unsigned int n_seg_min = 6u; n_seg_min >= displacedParameters.minNumberOfHits; --n_seg_min)
0122       lookForSegments(displacedParameters, n_seg_min, rechits, recHits_per_layer, used, segments);
0123   };
0124   // Not currently used

0125   // auto doWide = [&] () {

0126   //    for(unsigned int n_seg_min = 6u; n_seg_min >= wideParameters.minNumberOfHits; --n_seg_min)

0127   //        lookForSegments(wideParameters,n_seg_min,rechits,recHits_per_layer, used,segments);

0128   // };

0129   auto printSegments = [&] {
0130 #ifdef EDM_ML_DEBUG  // have lines below only compiled when in debug mode

0131     for (unsigned int iS = 0; iS < segments.size(); ++iS) {
0132       const auto& seg = segments[iS];
0133       ME0DetId chId(seg.me0DetId());
0134       const auto& rechits = seg.specificRecHits();
0135       edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegAlgoRU] segment in chamber " << chId << " which contains "
0136                                        << rechits.size() << " rechits and with specs: \n"
0137                                        << seg;
0138       for (auto rh = rechits.begin(); rh != rechits.end(); ++rh) {
0139         auto me0id = rh->me0Id();
0140         edm::LogVerbatim("ME0SegAlgoRU") << "[RecHit :: Loc x = " << std::showpos << std::setw(9)
0141                                          << rh->localPosition().x() << " Loc y = " << std::showpos << std::setw(9)
0142                                          << rh->localPosition().y() << " Time = " << std::showpos << rh->tof() << " -- "
0143                                          << me0id.rawId() << " = " << me0id << " ]";
0144       }
0145     }
0146 #endif
0147   };
0148 
0149   //If we arent doing collisions, do a single iteration

0150   if (!doCollisions) {
0151     doDisplaced();
0152     printSegments();
0153     return segments;
0154   }
0155 
0156   //Iteration 1: STD processing

0157   doStd();
0158 
0159   if (false) {
0160     //How the CSC algorithm ~does iterations. for now not considering

0161     //displaced muons will not worry about it  later

0162     //Iteration 2a: If we don't allow wide segments simply do displaced

0163     // Or if we already found a segment simply skip to displaced

0164     if (!allowWideSegments || !segments.empty()) {
0165       doDisplaced();
0166       return segments;
0167     }
0168     //doWide();

0169     doDisplaced();
0170   }
0171   printSegments();
0172   return segments;
0173 }
0174 
0175 void ME0SegAlgoRU::lookForSegments(const SegmentParameters& params,
0176                                    const unsigned int n_seg_min,
0177                                    const HitAndPositionContainer& rechits,
0178                                    const std::vector<unsigned int>& recHits_per_layer,
0179                                    BoolContainer& used,
0180                                    std::vector<ME0Segment>& segments) const {
0181   auto ib = rechits.begin();
0182   auto ie = rechits.end();
0183   std::vector<std::pair<float, HitAndPositionPtrContainer> > proto_segments;
0184   // the first hit is taken from the back

0185   for (auto i1 = ib; i1 != ie; ++i1) {
0186     const auto& h1 = *i1;
0187 
0188     //skip if rh is used and the layer tat has big rh multiplicity(>25RHs)

0189     if (used[h1.idx])
0190       continue;
0191     if (recHits_per_layer[h1.layer - 1] > 100)
0192       continue;
0193 
0194     //      bool segok = false;

0195     // the second hit from the front

0196     for (auto i2 = ie - 1; i2 != i1; --i2) {
0197       //            if(segok) break; //Currently turned off

0198       const auto& h2 = *i2;
0199 
0200       //skip if rh is used and the layer tat has big rh multiplicity(>25RHs)

0201       if (used[h2.idx])
0202         continue;
0203       if (recHits_per_layer[h2.layer - 1] > 100)
0204         continue;
0205 
0206       //Stop if the distance between layers is not large enough

0207       if ((std::abs(int(h2.layer) - int(h1.layer)) + 1) < int(n_seg_min))
0208         break;
0209 
0210       if (std::fabs(h1.rh->tof() - h2.rh->tof()) > params.maxTOFDiff + 1.0)
0211         continue;
0212       if (!areHitsCloseInEta(params.maxETASeeds, params.requireBeamConstr, h1.gp, h2.gp))
0213         continue;
0214       if (!areHitsCloseInGlobalPhi(params.maxPhiSeeds, abs(int(h2.layer) - int(h1.layer)), h1.gp, h2.gp))
0215         continue;
0216 
0217       HitAndPositionPtrContainer current_proto_segment;
0218       std::unique_ptr<MuonSegFit> current_fit;
0219       current_fit = addHit(current_proto_segment, h1);
0220       current_fit = addHit(current_proto_segment, h2);
0221 
0222       tryAddingHitsToSegment(params.maxTOFDiff,
0223                              params.maxETASeeds,
0224                              params.maxPhiAdditional,
0225                              params.maxChi2Additional,
0226                              current_fit,
0227                              current_proto_segment,
0228                              used,
0229                              i1,
0230                              i2);
0231 
0232       if (current_proto_segment.size() > n_seg_min)
0233         pruneBadHits(params.maxChi2Prune, current_proto_segment, current_fit, n_seg_min);
0234 
0235       edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegAlgoRU::lookForSegments] # of hits in segment "
0236                                        << current_proto_segment.size() << " min # " << n_seg_min << " => "
0237                                        << (current_proto_segment.size() >= n_seg_min) << " chi2/ndof "
0238                                        << current_fit->chi2() / current_fit->ndof() << " => "
0239                                        << (current_fit->chi2() / current_fit->ndof() < params.maxChi2GoodSeg)
0240                                        << std::endl;
0241 
0242       if (current_proto_segment.size() < n_seg_min)
0243         continue;
0244       const float current_metric = current_fit->chi2() / current_fit->ndof();
0245       if (current_metric > params.maxChi2GoodSeg)
0246         continue;
0247 
0248       if (params.requireCentralBX) {
0249         int nCentral = 0;
0250         int nNonCentral = 0;
0251         for (const auto* rh : current_proto_segment) {
0252           if (std::fabs(rh->rh->tof()) < 2)
0253             nCentral++;
0254           else
0255             nNonCentral++;
0256         }
0257         if (nNonCentral >= nCentral)
0258           continue;
0259       }
0260       //            segok = true;

0261 
0262       proto_segments.emplace_back(current_metric, current_proto_segment);
0263     }
0264   }
0265   addUniqueSegments(proto_segments, segments, used);
0266 }
0267 
0268 void ME0SegAlgoRU::addUniqueSegments(SegmentByMetricContainer& proto_segments,
0269                                      std::vector<ME0Segment>& segments,
0270                                      BoolContainer& used) const {
0271   std::sort(proto_segments.begin(),
0272             proto_segments.end(),
0273             [](const std::pair<float, HitAndPositionPtrContainer>& a,
0274                const std::pair<float, HitAndPositionPtrContainer>& b) { return a.first < b.first; });
0275 
0276   //Now add to the collect based on minChi2 marking the hits as used after

0277   std::vector<unsigned int> usedHits;
0278   for (unsigned int iS = 0; iS < proto_segments.size(); ++iS) {
0279     HitAndPositionPtrContainer currentProtoSegment = proto_segments[iS].second;
0280 
0281     //check to see if we already used thes hits this round

0282     bool alreadyFilled = false;
0283     for (unsigned int iH = 0; iH < currentProtoSegment.size(); ++iH) {
0284       for (unsigned int iOH = 0; iOH < usedHits.size(); ++iOH) {
0285         if (usedHits[iOH] != currentProtoSegment[iH]->idx)
0286           continue;
0287         alreadyFilled = true;
0288         break;
0289       }
0290     }
0291     if (alreadyFilled)
0292       continue;
0293     for (const auto* h : currentProtoSegment) {
0294       usedHits.push_back(h->idx);
0295       used[h->idx] = true;
0296     }
0297 
0298     std::unique_ptr<MuonSegFit> current_fit = makeFit(currentProtoSegment);
0299 
0300     // Create an actual ME0Segment - retrieve all info from the fit

0301     // calculate the timing fron rec hits associated to the TrackingRecHits used

0302     // to fit the segment

0303     float averageTime = 0.;
0304     for (const auto* h : currentProtoSegment) {
0305       averageTime += h->rh->tof();
0306     }
0307     averageTime /= float(currentProtoSegment.size());
0308     float timeUncrt = 0.;
0309     for (const auto* h : currentProtoSegment) {
0310       timeUncrt += (h->rh->tof() - averageTime) * (h->rh->tof() - averageTime);
0311     }
0312     timeUncrt = std::sqrt(timeUncrt / float(currentProtoSegment.size() - 1));
0313 
0314     std::sort(currentProtoSegment.begin(),
0315               currentProtoSegment.end(),
0316               [](const HitAndPosition* a, const HitAndPosition* b) { return a->layer < b->layer; });
0317 
0318     std::vector<const ME0RecHit*> bareRHs;
0319     bareRHs.reserve(currentProtoSegment.size());
0320     for (const auto* rh : currentProtoSegment)
0321       bareRHs.push_back(rh->rh);
0322     const float dPhi = theChamber->computeDeltaPhi(current_fit->intercept(), current_fit->localdir());
0323     ME0Segment temp(bareRHs,
0324                     current_fit->intercept(),
0325                     current_fit->localdir(),
0326                     current_fit->covarianceMatrix(),
0327                     current_fit->chi2(),
0328                     averageTime,
0329                     timeUncrt,
0330                     dPhi);
0331     segments.push_back(temp);
0332   }
0333 }
0334 
0335 void ME0SegAlgoRU::tryAddingHitsToSegment(const float maxTOF,
0336                                           const float maxETA,
0337                                           const float maxPhi,
0338                                           const float maxChi2,
0339                                           std::unique_ptr<MuonSegFit>& current_fit,
0340                                           HitAndPositionPtrContainer& proto_segment,
0341                                           const BoolContainer& used,
0342                                           HitAndPositionContainer::const_iterator i1,
0343                                           HitAndPositionContainer::const_iterator i2) const {
0344   // Iterate over the layers with hits in the chamber

0345   // Skip the layers containing the segment endpoints

0346   // Test each hit on the other layers to see if it is near the segment

0347   // If it is, see whether there is already a hit on the segment from the same layer

0348   //    - if so, and there are more than 2 hits on the segment, copy the segment,

0349   //      replace the old hit with the new hit. If the new segment chi2 is better

0350   //      then replace the original segment with the new one (by swap)

0351   //    - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory

0352   //      then replace the original segment with the new one (by swap)

0353 
0354   //Hits are ordered by layer, "i1" is the inner hit and i2 is the outer hit

0355   //so possible hits to add must be between these two iterators

0356   for (auto iH = i1 + 1; iH != i2; ++iH) {
0357     if (iH->layer == i1->layer)
0358       continue;
0359     if (iH->layer == i2->layer)
0360       break;
0361     if (used[iH->idx])
0362       continue;
0363     if (!areHitsConsistentInTime(maxTOF, proto_segment, *iH))
0364       continue;
0365     if (!isHitNearSegment(maxETA, maxPhi, current_fit, proto_segment, *iH))
0366       continue;
0367     if (hasHitOnLayer(proto_segment, iH->layer))
0368       compareProtoSegment(current_fit, proto_segment, *iH);
0369     else
0370       increaseProtoSegment(maxChi2, current_fit, proto_segment, *iH);
0371   }
0372 }
0373 
0374 bool ME0SegAlgoRU::areHitsCloseInEta(const float maxETA,
0375                                      const bool beamConst,
0376                                      const GlobalPoint& h1,
0377                                      const GlobalPoint& h2) const {
0378   float diff = std::fabs(h1.eta() - h2.eta());
0379   edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegAlgoRU::areHitsCloseInEta] gp1 = " << h1 << " in eta part = " << h1.eta()
0380                                    << " and gp2 = " << h2 << " in eta part = " << h2.eta() << " ==> dEta = " << diff
0381                                    << " ==> return " << (diff < 0.1) << std::endl;
0382   //temp for floating point comparision...maxEta is the difference between partitions, so x1.5 to take into account non-circle geom.

0383   return (diff < std::max(maxETA, 0.01f));
0384 }
0385 
0386 bool ME0SegAlgoRU::areHitsCloseInGlobalPhi(const float maxPHI,
0387                                            const unsigned int nLayDisp,
0388                                            const GlobalPoint& h1,
0389                                            const GlobalPoint& h2) const {
0390   float dphi12 = deltaPhi(h1.barePhi(), h2.barePhi());
0391   edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegAlgoRU::areHitsCloseInGlobalPhi] gp1 = " << h1 << " and gp2 = " << h2
0392                                    << " ==> dPhi = " << dphi12 << " ==> return "
0393                                    << (fabs(dphi12) < std::max(maxPHI, 0.02f)) << std::endl;
0394   return fabs(dphi12) < std::max(maxPHI, float(float(nLayDisp) * 0.004));
0395 }
0396 
0397 bool ME0SegAlgoRU::isHitNearSegment(const float maxETA,
0398                                     const float maxPHI,
0399                                     const std::unique_ptr<MuonSegFit>& fit,
0400                                     const HitAndPositionPtrContainer& proto_segment,
0401                                     const HitAndPosition& h) const {
0402   //Get average eta, based on the two seeds...asssumes that we have not started pruning yet!

0403   const float avgETA = (proto_segment[1]->gp.eta() + proto_segment[0]->gp.eta()) / 2.;
0404   //    edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegAlgoRU::isHitNearSegment] average eta = "<<avgETA<<" additionalHit eta = "<<h.gp.eta() <<" ==> dEta = "<<std::fabs(h.gp.eta() - avgETA) <<" ==> return "<<(std::fabs(h.gp.eta() - avgETA) < std::max(maxETA,0.01f ))<<std::endl;

0405   if (std::fabs(h.gp.eta() - avgETA) > std::max(maxETA, 0.01f))
0406     return false;
0407 
0408   //Now check the dPhi based on the segment fit

0409   GlobalPoint globIntercept = globalAtZ(fit, h.lp.z());
0410   float dPhi = deltaPhi(h.gp.barePhi(), globIntercept.phi());
0411   //check to see if it is inbetween the two rolls of the outer and inner hits

0412   //    edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegAlgoRU::isHitNearSegment] projected phi = "<<globIntercept.phi()<<" additionalHit phi = "<<h.gp.barePhi() <<" ==> dPhi = "<<dPhi <<" ==> return "<<(std::fabs(dPhi) <  std::max(maxPHI,0.001f))<<std::endl;

0413   return (std::fabs(dPhi) < std::max(maxPHI, 0.001f));
0414 }
0415 
0416 bool ME0SegAlgoRU::areHitsConsistentInTime(const float maxTOF,
0417                                            const HitAndPositionPtrContainer& proto_segment,
0418                                            const HitAndPosition& h) const {
0419   std::vector<float> tofs;
0420   tofs.reserve(proto_segment.size() + 1);
0421   tofs.push_back(h.rh->tof());
0422   for (const auto* ih : proto_segment)
0423     tofs.push_back(ih->rh->tof());
0424   std::sort(tofs.begin(), tofs.end());
0425   if (std::fabs(tofs.front() - tofs.back()) < maxTOF + 1.0)
0426     return true;
0427   return false;
0428 }
0429 
0430 GlobalPoint ME0SegAlgoRU::globalAtZ(const std::unique_ptr<MuonSegFit>& fit, float z) const {
0431   float x = fit->xfit(z);
0432   float y = fit->yfit(z);
0433   return theChamber->toGlobal(LocalPoint(x, y, z));
0434 }
0435 
0436 std::unique_ptr<MuonSegFit> ME0SegAlgoRU::addHit(HitAndPositionPtrContainer& proto_segment,
0437                                                  const HitAndPosition& aHit) const {
0438   proto_segment.push_back(&aHit);
0439   // make a fit

0440   return makeFit(proto_segment);
0441 }
0442 
0443 std::unique_ptr<MuonSegFit> ME0SegAlgoRU::makeFit(const HitAndPositionPtrContainer& proto_segment) const {
0444   // for ME0 we take the me0rechit from the proto_segment we transform into Tracking Rechits

0445   // the local rest frame is the ME0Chamber

0446   MuonSegFit::MuonRecHitContainer muonRecHits;
0447   for (auto rh = proto_segment.begin(); rh < proto_segment.end(); rh++) {
0448     ME0RecHit* newRH = (*rh)->rh->clone();
0449     newRH->setPosition((*rh)->lp);
0450     MuonSegFit::MuonRecHitPtr trkRecHit(newRH);
0451     muonRecHits.push_back(trkRecHit);
0452   }
0453   std::unique_ptr<MuonSegFit> currentFit(new MuonSegFit(muonRecHits));
0454   currentFit->fit();
0455   return currentFit;
0456 }
0457 
0458 void ME0SegAlgoRU::pruneBadHits(const float maxChi2,
0459                                 HitAndPositionPtrContainer& proto_segment,
0460                                 std::unique_ptr<MuonSegFit>& fit,
0461                                 const unsigned int n_seg_min) const {
0462   while (proto_segment.size() > n_seg_min && fit->chi2() / fit->ndof() > maxChi2) {
0463     float maxDev = -1;
0464     HitAndPositionPtrContainer::iterator worstHit;
0465     for (auto it = proto_segment.begin(); it != proto_segment.end(); ++it) {
0466       const float dev = getHitSegChi2(fit, *(*it)->rh);
0467       if (dev < maxDev)
0468         continue;
0469       maxDev = dev;
0470       worstHit = it;
0471     }
0472     edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegAlgoRU::pruneBadHits] pruning one hit-> layer: " << (*worstHit)->layer
0473                                      << " eta: " << (*worstHit)->gp.eta() << " phi: " << (*worstHit)->gp.phi()
0474                                      << " old chi2/dof: " << fit->chi2() / fit->ndof() << std::endl;
0475     proto_segment.erase(worstHit);
0476     fit = makeFit(proto_segment);
0477   }
0478 }
0479 
0480 float ME0SegAlgoRU::getHitSegChi2(const std::unique_ptr<MuonSegFit>& fit, const ME0RecHit& hit) const {
0481   const auto lp = hit.localPosition();
0482   const auto le = hit.localPositionError();
0483   const float du = fit->xdev(lp.x(), lp.z());
0484   const float dv = fit->ydev(lp.y(), lp.z());
0485 
0486   ROOT::Math::SMatrix<double, 2, 2, ROOT::Math::MatRepSym<double, 2> > IC;
0487   IC(0, 0) = le.xx();
0488   IC(1, 0) = le.xy();
0489   IC(1, 1) = le.yy();
0490 
0491   // Invert covariance matrix

0492   IC.Invert();
0493   return du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
0494 }
0495 
0496 bool ME0SegAlgoRU::hasHitOnLayer(const HitAndPositionPtrContainer& proto_segment, const unsigned int layer) const {
0497   for (const auto* h : proto_segment)
0498     if (h->layer == layer)
0499       return true;
0500   return false;
0501 }
0502 
0503 void ME0SegAlgoRU::compareProtoSegment(std::unique_ptr<MuonSegFit>& current_fit,
0504                                        HitAndPositionPtrContainer& current_proto_segment,
0505                                        const HitAndPosition& new_hit) const {
0506   const HitAndPosition* old_hit = nullptr;
0507   HitAndPositionPtrContainer new_proto_segment = current_proto_segment;
0508 
0509   for (auto it = new_proto_segment.begin(); it != new_proto_segment.end();) {
0510     if ((*it)->layer == new_hit.layer) {
0511       old_hit = *it;
0512       it = new_proto_segment.erase(it);
0513     } else {
0514       ++it;
0515     }
0516   }
0517   if (old_hit == nullptr)
0518     return;
0519   auto new_fit = addHit(new_proto_segment, new_hit);
0520 
0521   //If on the same strip but different BX choose the closest

0522   bool useNew = false;
0523   if (old_hit->lp == new_hit.lp) {
0524     float avgtof = 0;
0525     for (const auto* h : current_proto_segment)
0526       if (old_hit != h)
0527         avgtof += h->rh->tof();
0528     avgtof /= float(current_proto_segment.size() - 1);
0529     if (std::abs(avgtof - new_hit.rh->tof()) < std::abs(avgtof - old_hit->rh->tof()))
0530       useNew = true;
0531   }  //otherwise base it on chi2

0532   else if (new_fit->chi2() < current_fit->chi2())
0533     useNew = true;
0534 
0535   if (useNew) {
0536     current_proto_segment = new_proto_segment;
0537     current_fit = std::move(new_fit);
0538   }
0539 }
0540 
0541 void ME0SegAlgoRU::increaseProtoSegment(const float maxChi2,
0542                                         std::unique_ptr<MuonSegFit>& current_fit,
0543                                         HitAndPositionPtrContainer& current_proto_segment,
0544                                         const HitAndPosition& new_hit) const {
0545   HitAndPositionPtrContainer new_proto_segment = current_proto_segment;
0546   auto new_fit = addHit(new_proto_segment, new_hit);
0547   //    edm::LogVerbatim("ME0SegAlgoRU") << "[ME0SegAlgoRU::increaseProtoSegment] new chi2 = "<<new_fit->chi2()<<" new chi2/dof = "<<new_fit->chi2()/new_fit->ndof() <<" ==> add: " << (new_fit->chi2()/new_fit->ndof() < maxChi2 ) <<std::endl;

0548   if (new_fit->chi2() / new_fit->ndof() < maxChi2) {
0549     current_proto_segment = new_proto_segment;
0550     current_fit = std::move(new_fit);
0551   }
0552 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!