Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 /**
0002  * \file GE0SegAlgRU.cc
0003  *
0004  *  \author I. J. Watson adaptations for GE0
0005  *  \author M. Maggi for ME0
0006  *  \from V.Palichik & N.Voytishin
0007  *  \some functions and structure taken from SK algo by M.Sani and SegFit class by T.Cox
0008  */
0009 #include "GE0SegAlgoRU.h"
0010 #include "MuonSegFit.h"
0011 #include "Geometry/GEMGeometry/interface/GEMChamber.h"
0012 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0013 #include "DataFormats/Math/interface/deltaPhi.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016 
0017 #include <algorithm>
0018 #include <cmath>
0019 #include <iostream>
0020 #include <string>
0021 
0022 #include <Math/Functions.h>
0023 #include <Math/SVector.h>
0024 #include <Math/SMatrix.h>
0025 
0026 GE0SegAlgoRU::GE0SegAlgoRU(const edm::ParameterSet& ps) : GEMSegmentAlgorithmBase(ps), myName("GE0SegAlgoRU") {
0027   allowWideSegments = ps.getParameter<bool>("allowWideSegments");
0028   doCollisions = ps.getParameter<bool>("doCollisions");
0029 
0030   stdParameters.maxChi2Additional = ps.getParameter<double>("maxChi2Additional");
0031   stdParameters.maxChi2Prune = ps.getParameter<double>("maxChi2Prune");
0032   stdParameters.maxChi2GoodSeg = ps.getParameter<double>("maxChi2GoodSeg");
0033   stdParameters.maxPhiSeeds = ps.getParameter<double>("maxPhiSeeds");
0034   stdParameters.maxPhiAdditional = ps.getParameter<double>("maxPhiAdditional");
0035   stdParameters.maxETASeeds = ps.getParameter<double>("maxETASeeds");
0036   stdParameters.requireCentralBX = ps.getParameter<bool>("requireCentralBX");
0037   stdParameters.minNumberOfHits = ps.getParameter<unsigned int>("minNumberOfHits");
0038   stdParameters.maxNumberOfHits = ps.getParameter<unsigned int>("maxNumberOfHits");
0039   stdParameters.maxNumberOfHitsPerLayer = ps.getParameter<unsigned int>("maxNumberOfHitsPerLayer");
0040   stdParameters.requireBeamConstr = true;
0041 
0042   wideParameters = stdParameters;
0043   wideParameters.maxChi2Prune *= 2;
0044   wideParameters.maxChi2GoodSeg *= 2;
0045   wideParameters.maxPhiSeeds *= 2;
0046   wideParameters.maxPhiAdditional *= 2;
0047   wideParameters.minNumberOfHits += 1;
0048 
0049   displacedParameters = stdParameters;
0050   displacedParameters.maxChi2Additional *= 2;
0051   displacedParameters.maxChi2Prune = 100;
0052   displacedParameters.maxChi2GoodSeg *= 5;
0053   displacedParameters.maxPhiSeeds *= 2;
0054   displacedParameters.maxPhiAdditional *= 2;
0055   displacedParameters.maxETASeeds *= 2;
0056   displacedParameters.requireBeamConstr = false;
0057 
0058   LogDebug("GE0SegAlgoRU") << myName << " has algorithm cuts set to: \n"
0059                            << "--------------------------------------------------------------------\n"
0060                            << "allowWideSegments   = " << allowWideSegments << "\n"
0061                            << "doCollisions        = " << doCollisions << "\n"
0062                            << "maxChi2Additional   = " << stdParameters.maxChi2Additional << "\n"
0063                            << "maxChi2Prune        = " << stdParameters.maxChi2Prune << "\n"
0064                            << "maxChi2GoodSeg      = " << stdParameters.maxChi2GoodSeg << "\n"
0065                            << "maxPhiSeeds         = " << stdParameters.maxPhiSeeds << "\n"
0066                            << "maxPhiAdditional    = " << stdParameters.maxPhiAdditional << "\n"
0067                            << "maxETASeeds         = " << stdParameters.maxETASeeds << "\n"
0068                            << std::endl;
0069 
0070   theChamber = nullptr;
0071 }
0072 
0073 std::vector<GEMSegment> GE0SegAlgoRU::run(const GEMEnsemble& ensemble, const std::vector<const GEMRecHit*>& rechits) {
0074   HitAndPositionContainer hitAndPositions;
0075   auto& superchamber = ensemble.first;
0076   for (const auto& rechit : rechits) {
0077     const GEMEtaPartition* part = ensemble.second.at(rechit->gemId().rawId());
0078     GlobalPoint glb = part->toGlobal(rechit->localPosition());
0079     LocalPoint nLoc = superchamber->toLocal(glb);
0080     hitAndPositions.emplace_back(&(*rechit), nLoc, glb, hitAndPositions.size());
0081   }
0082 
0083   LogDebug("GE0Segment|GE0") << "found " << hitAndPositions.size() << " rechits in superchamber " << superchamber->id();
0084   //sort by layer
0085   float z1 = superchamber->chamber(1)->position().z();
0086   float zlast = superchamber->chamber(superchamber->nChambers())->position().z();
0087   if (z1 < zlast)
0088     std::sort(hitAndPositions.begin(), hitAndPositions.end(), [](const HitAndPosition& h1, const HitAndPosition& h2) {
0089       return h1.layer < h2.layer;
0090     });
0091   else
0092     std::sort(hitAndPositions.begin(), hitAndPositions.end(), [](const HitAndPosition& h1, const HitAndPosition& h2) {
0093       return h1.layer > h2.layer;
0094     });
0095   return run(ensemble.first, hitAndPositions);
0096 }
0097 
0098 std::vector<GEMSegment> GE0SegAlgoRU::run(const GEMSuperChamber* chamber, const HitAndPositionContainer& rechits) {
0099 #ifdef EDM_ML_DEBUG  // have lines below only compiled when in debug mode
0100   GEMDetId chId(chamber->id());
0101   edm::LogVerbatim("GE0SegAlgoRU") << "[GEMSegmentAlgorithm::run] build segments in chamber " << chId
0102                                    << " which contains " << rechits.size() << " rechits";
0103   for (const auto& h : rechits) {
0104     auto ge0id = h.rh->gemId();
0105     auto rhLP = h.lp;
0106     edm::LogVerbatim("GE0SegAlgoRU") << "[RecHit :: Loc x = " << std::showpos << std::setw(9) << rhLP.x()
0107                                      << " Glb y = " << std::showpos << std::setw(9) << rhLP.y()
0108                                      << " Time = " << std::showpos << h.rh->BunchX() << " -- " << ge0id.rawId() << " = "
0109                                      << ge0id << " ]" << std::endl;
0110   }
0111 #endif
0112 
0113   if (rechits.size() < stdParameters.minNumberOfHits || rechits.size() > stdParameters.maxNumberOfHits) {
0114     return std::vector<GEMSegment>();
0115   }
0116 
0117   theChamber = chamber;
0118 
0119   std::vector<unsigned int> recHits_per_layer(theChamber->nChambers(), 0);
0120   for (const auto& rechit : rechits) {
0121     recHits_per_layer[rechit.layer - 1]++;
0122   }
0123 
0124   BoolContainer used(rechits.size(), false);
0125 
0126   // We have at least 2 hits. We intend to try all possible pairs of hits to start
0127   // segment building. 'All possible' means each hit lies on different layers in the chamber.
0128   // after all same size segs are build we get rid of the overcrossed segments using the chi2 criteria
0129   // the hits from the segs that are left are marked as used and are not added to segs in future iterations
0130   // the hits from 3p segs are marked as used separately in order to try to assamble them in longer segments
0131   // in case there is a second pass
0132 
0133   // Choose first hit (as close to IP as possible) h1 and second hit
0134   // (as far from IP as possible) h2 To do this we iterate over hits
0135   // in the chamber by layer - pick two layers.  Then we
0136   // iterate over hits within each of these layers and pick h1 and h2
0137   // these.  If they are 'close enough' we build an empty
0138   // segment.  Then try adding hits to this segment.
0139 
0140   std::vector<GEMSegment> segments;
0141 
0142   auto doStd = [&]() {
0143     for (unsigned int n_seg_min = 6u; n_seg_min >= stdParameters.minNumberOfHits; --n_seg_min)
0144       lookForSegments(stdParameters, n_seg_min, rechits, recHits_per_layer, used, segments);
0145   };
0146   auto doDisplaced = [&]() {
0147     for (unsigned int n_seg_min = 6u; n_seg_min >= displacedParameters.minNumberOfHits; --n_seg_min)
0148       lookForSegments(displacedParameters, n_seg_min, rechits, recHits_per_layer, used, segments);
0149   };
0150   // Not currently used
0151   // auto doWide = [&] () {
0152   //    for(unsigned int n_seg_min = 6u; n_seg_min >= wideParameters.minNumberOfHits; --n_seg_min)
0153   //        lookForSegments(wideParameters,n_seg_min,rechits,recHits_per_layer, used,segments);
0154   // };
0155   auto printSegments = [&] {
0156 #ifdef EDM_ML_DEBUG  // have lines below only compiled when in debug mode
0157     for (const auto& seg : segments) {
0158       GEMDetId chId(seg.gemDetId());
0159       const auto& rechits = seg.specificRecHits();
0160       edm::LogVerbatim("GE0SegAlgoRU") << "[GE0SegAlgoRU] segment in chamber " << chId << " which contains "
0161                                        << rechits.size() << " rechits and with specs: \n"
0162                                        << seg;
0163       for (const auto& rh : rechits) {
0164         auto ge0id = rh.gemId();
0165         edm::LogVerbatim("GE0SegAlgoRU") << "[RecHit :: Loc x = " << std::showpos << std::setw(9)
0166                                          << rh.localPosition().x() << " Loc y = " << std::showpos << std::setw(9)
0167                                          << rh.localPosition().y() << " Time = " << std::showpos << rh.BunchX()
0168                                          << " -- " << ge0id.rawId() << " = " << ge0id << " ]";
0169       }
0170     }
0171 #endif
0172   };
0173 
0174   //If we arent doing collisions, do a single iteration
0175   if (!doCollisions) {
0176     doDisplaced();
0177     printSegments();
0178     return segments;
0179   }
0180 
0181   //Iteration 1: STD processing
0182   doStd();
0183 
0184   if (false) {
0185     //How the CSC algorithm ~does iterations. for now not considering
0186     //displaced muons will not worry about it  later
0187     //Iteration 2a: If we don't allow wide segments simply do displaced
0188     // Or if we already found a segment simply skip to displaced
0189     if (!allowWideSegments || !segments.empty()) {
0190       doDisplaced();
0191       return segments;
0192     }
0193     //doWide();
0194     doDisplaced();
0195   }
0196   printSegments();
0197   return segments;
0198 }
0199 
0200 void GE0SegAlgoRU::lookForSegments(const SegmentParameters& params,
0201                                    const unsigned int n_seg_min,
0202                                    const HitAndPositionContainer& rechits,
0203                                    const std::vector<unsigned int>& recHits_per_layer,
0204                                    BoolContainer& used,
0205                                    std::vector<GEMSegment>& segments) const {
0206   auto ib = rechits.begin();
0207   auto ie = rechits.end();
0208   std::vector<std::pair<float, HitAndPositionPtrContainer> > proto_segments;
0209   // the first hit is taken from the back
0210   for (auto i1 = ib; i1 != ie; ++i1) {
0211     const auto& h1 = *i1;
0212 
0213     //skip if rh is used and the layer tat has big rh multiplicity(>25RHs)
0214     if (used[h1.idx])
0215       continue;
0216     if (recHits_per_layer[h1.layer - 1] > params.maxNumberOfHitsPerLayer)
0217       continue;
0218 
0219     // the second hit from the front
0220     for (auto i2 = ie - 1; i2 != i1; --i2) {
0221       const auto& h2 = *i2;
0222 
0223       //skip if rh is used and the layer tat has big rh multiplicity(>25RHs)
0224       if (used[h2.idx])
0225         continue;
0226       if (recHits_per_layer[h2.layer - 1] > params.maxNumberOfHitsPerLayer)
0227         continue;
0228 
0229       //Stop if the distance between layers is not large enough
0230       if ((std::abs(int(h2.layer) - int(h1.layer)) + 1) < int(n_seg_min))
0231         break;
0232 
0233       if (!areHitsCloseInEta(params.maxETASeeds, params.requireBeamConstr, h1.gp, h2.gp))
0234         continue;
0235       if (!areHitsCloseInGlobalPhi(params.maxPhiSeeds, std::abs(int(h2.layer) - int(h1.layer)), h1.gp, h2.gp))
0236         continue;
0237 
0238       HitAndPositionPtrContainer current_proto_segment;
0239       std::unique_ptr<MuonSegFit> current_fit;
0240       current_fit = addHit(current_proto_segment, h1);
0241       current_fit = addHit(current_proto_segment, h2);
0242 
0243       tryAddingHitsToSegment(params.maxETASeeds,
0244                              params.maxPhiAdditional,
0245                              params.maxChi2Additional,
0246                              current_fit,
0247                              current_proto_segment,
0248                              used,
0249                              i1,
0250                              i2);
0251 
0252       if (current_proto_segment.size() > n_seg_min)
0253         pruneBadHits(params.maxChi2Prune, current_proto_segment, current_fit, n_seg_min);
0254 
0255       LogDebug("GE0SegAlgoRU") << "[GE0SegAlgoRU::lookForSegments] # of hits in segment "
0256                                << current_proto_segment.size() << " min # " << n_seg_min << " => "
0257                                << (current_proto_segment.size() >= n_seg_min) << " chi2/ndof "
0258                                << current_fit->chi2() / current_fit->ndof() << " => "
0259                                << (current_fit->chi2() / current_fit->ndof() < params.maxChi2GoodSeg) << std::endl;
0260 
0261       if (current_proto_segment.size() < n_seg_min)
0262         continue;
0263       const float current_metric = current_fit->chi2() / current_fit->ndof();
0264       if (current_metric > params.maxChi2GoodSeg)
0265         continue;
0266 
0267       if (params.requireCentralBX) {
0268         int nCentral = 0;
0269         int nNonCentral = 0;
0270         for (const auto* rh : current_proto_segment) {
0271           if (std::abs(rh->rh->BunchX()) < 2)
0272             nCentral++;
0273           else
0274             nNonCentral++;
0275         }
0276         if (nNonCentral >= nCentral)
0277           continue;
0278       }
0279 
0280       proto_segments.emplace_back(current_metric, current_proto_segment);
0281     }
0282   }
0283   addUniqueSegments(proto_segments, segments, used);
0284 }
0285 
0286 void GE0SegAlgoRU::addUniqueSegments(SegmentByMetricContainer& proto_segments,
0287                                      std::vector<GEMSegment>& segments,
0288                                      BoolContainer& used) const {
0289   std::sort(proto_segments.begin(),
0290             proto_segments.end(),
0291             [](const std::pair<float, HitAndPositionPtrContainer>& a,
0292                const std::pair<float, HitAndPositionPtrContainer>& b) { return a.first < b.first; });
0293 
0294   //Now add to the collect based on minChi2 marking the hits as used after
0295   std::vector<unsigned int> usedHits;
0296   for (auto& container : proto_segments) {
0297     HitAndPositionPtrContainer currentProtoSegment = container.second;
0298 
0299     //check to see if we already used thes hits this round
0300     bool alreadyFilled = false;
0301     for (const auto& h : currentProtoSegment) {
0302       for (unsigned int iOH = 0; iOH < usedHits.size(); ++iOH) {
0303         if (usedHits[iOH] != h->idx)
0304           continue;
0305         alreadyFilled = true;
0306         break;
0307       }
0308     }
0309     if (alreadyFilled)
0310       continue;
0311     for (const auto* h : currentProtoSegment) {
0312       usedHits.push_back(h->idx);
0313       used[h->idx] = true;
0314     }
0315 
0316     std::unique_ptr<MuonSegFit> current_fit = makeFit(currentProtoSegment);
0317 
0318     // Create an actual GEMSegment - retrieve all info from the fit
0319     // calculate the timing fron rec hits associated to the TrackingRecHits used
0320     // to fit the segment
0321     float averageBX = 0.;
0322     for (const auto* h : currentProtoSegment) {
0323       averageBX += h->rh->BunchX();
0324     }
0325     averageBX /= int(currentProtoSegment.size());
0326 
0327     std::sort(currentProtoSegment.begin(),
0328               currentProtoSegment.end(),
0329               [](const HitAndPosition* a, const HitAndPosition* b) { return a->layer < b->layer; });
0330 
0331     std::vector<const GEMRecHit*> bareRHs;
0332     bareRHs.reserve(currentProtoSegment.size());
0333     for (const auto* rh : currentProtoSegment)
0334       bareRHs.push_back(rh->rh);
0335     const float dPhi = theChamber->computeDeltaPhi(current_fit->intercept(), current_fit->localdir());
0336     GEMSegment temp(bareRHs,
0337                     current_fit->intercept(),
0338                     current_fit->localdir(),
0339                     current_fit->covarianceMatrix(),
0340                     current_fit->chi2(),
0341                     averageBX,
0342                     dPhi);
0343     segments.push_back(temp);
0344   }
0345 }
0346 
0347 void GE0SegAlgoRU::tryAddingHitsToSegment(const float maxETA,
0348                                           const float maxPhi,
0349                                           const float maxChi2,
0350                                           std::unique_ptr<MuonSegFit>& current_fit,
0351                                           HitAndPositionPtrContainer& proto_segment,
0352                                           const BoolContainer& used,
0353                                           HitAndPositionContainer::const_iterator i1,
0354                                           HitAndPositionContainer::const_iterator i2) const {
0355   // Iterate over the layers with hits in the chamber
0356   // Skip the layers containing the segment endpoints
0357   // Test each hit on the other layers to see if it is near the segment
0358   // If it is, see whether there is already a hit on the segment from the same layer
0359   //    - if so, and there are more than 2 hits on the segment, copy the segment,
0360   //      replace the old hit with the new hit. If the new segment chi2 is better
0361   //      then replace the original segment with the new one (by swap)
0362   //    - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory
0363   //      then replace the original segment with the new one (by swap)
0364 
0365   //Hits are ordered by layer, "i1" is the inner hit and i2 is the outer hit
0366   //so possible hits to add must be between these two iterators
0367   for (auto iH = i1 + 1; iH != i2; ++iH) {
0368     if (iH->layer == i1->layer)
0369       continue;
0370     if (iH->layer == i2->layer)
0371       break;
0372     if (used[iH->idx])
0373       continue;
0374     if (!isHitNearSegment(maxETA, maxPhi, current_fit, proto_segment, *iH))
0375       continue;
0376     if (hasHitOnLayer(proto_segment, iH->layer))
0377       compareProtoSegment(current_fit, proto_segment, *iH);
0378     else
0379       increaseProtoSegment(maxChi2, current_fit, proto_segment, *iH);
0380   }
0381 }
0382 
0383 bool GE0SegAlgoRU::areHitsCloseInEta(const float maxETA,
0384                                      const bool beamConst,
0385                                      const GlobalPoint& h1,
0386                                      const GlobalPoint& h2) const {
0387   float diff = std::abs(h1.eta() - h2.eta());
0388   LogDebug("GE0SegAlgoRU") << "[GE0SegAlgoRU::areHitsCloseInEta] gp1 = " << h1 << " in eta part = " << h1.eta()
0389                            << " and gp2 = " << h2 << " in eta part = " << h2.eta() << " ==> dEta = " << diff
0390                            << " ==> return " << (diff < 0.1) << std::endl;
0391   //temp for floating point comparision...maxEta is the difference between partitions, so x1.5 to take into account non-circle geom.
0392   return (diff < std::max(maxETA, 0.01f));
0393 }
0394 
0395 bool GE0SegAlgoRU::areHitsCloseInGlobalPhi(const float maxPHI,
0396                                            const unsigned int nLayDisp,
0397                                            const GlobalPoint& h1,
0398                                            const GlobalPoint& h2) const {
0399   float dphi12 = deltaPhi(h1.barePhi(), h2.barePhi());
0400   LogDebug("GE0SegAlgoRU") << "[GE0SegAlgoRU::areHitsCloseInGlobalPhi] gp1 = " << h1 << " and gp2 = " << h2
0401                            << " ==> dPhi = " << dphi12 << " ==> return " << (std::abs(dphi12) < std::max(maxPHI, 0.02f))
0402                            << std::endl;
0403   return std::abs(dphi12) < std::max(maxPHI, float(float(nLayDisp) * 0.004));
0404 }
0405 
0406 bool GE0SegAlgoRU::isHitNearSegment(const float maxETA,
0407                                     const float maxPHI,
0408                                     const std::unique_ptr<MuonSegFit>& fit,
0409                                     const HitAndPositionPtrContainer& proto_segment,
0410                                     const HitAndPosition& h) const {
0411   //Get average eta, based on the two seeds...asssumes that we have not started pruning yet!
0412   const float avgETA = (proto_segment[1]->gp.eta() + proto_segment[0]->gp.eta()) / 2.;
0413   if (std::abs(h.gp.eta() - avgETA) > std::max(maxETA, 0.01f))
0414     return false;
0415 
0416   //Now check the dPhi based on the segment fit
0417   GlobalPoint globIntercept = globalAtZ(fit, h.lp.z());
0418   float dPhi = deltaPhi(h.gp.barePhi(), globIntercept.phi());
0419   //check to see if it is inbetween the two rolls of the outer and inner hits
0420   return (std::abs(dPhi) < std::max(maxPHI, 0.001f));
0421 }
0422 
0423 GlobalPoint GE0SegAlgoRU::globalAtZ(const std::unique_ptr<MuonSegFit>& fit, float z) const {
0424   float x = fit->xfit(z);
0425   float y = fit->yfit(z);
0426   return theChamber->toGlobal(LocalPoint(x, y, z));
0427 }
0428 
0429 std::unique_ptr<MuonSegFit> GE0SegAlgoRU::addHit(HitAndPositionPtrContainer& proto_segment,
0430                                                  const HitAndPosition& aHit) const {
0431   proto_segment.push_back(&aHit);
0432   // make a fit
0433   return makeFit(proto_segment);
0434 }
0435 
0436 std::unique_ptr<MuonSegFit> GE0SegAlgoRU::makeFit(const HitAndPositionPtrContainer& proto_segment) const {
0437   // for GE0 we take the gemrechit from the proto_segment we transform into Tracking Rechits
0438   // the local rest frame is the GEMSuperChamber
0439   MuonSegFit::MuonRecHitContainer muonRecHits;
0440   for (const auto& rh : proto_segment) {
0441     GEMRecHit* newRH = rh->rh->clone();
0442     newRH->setPosition(rh->lp);
0443     MuonSegFit::MuonRecHitPtr trkRecHit(newRH);
0444     muonRecHits.push_back(trkRecHit);
0445   }
0446   auto currentFit = std::make_unique<MuonSegFit>(muonRecHits);
0447   currentFit->fit();
0448   return currentFit;
0449 }
0450 
0451 void GE0SegAlgoRU::pruneBadHits(const float maxChi2,
0452                                 HitAndPositionPtrContainer& proto_segment,
0453                                 std::unique_ptr<MuonSegFit>& fit,
0454                                 const unsigned int n_seg_min) const {
0455   while (proto_segment.size() > n_seg_min && fit->chi2() / fit->ndof() > maxChi2) {
0456     float maxDev = -1;
0457     HitAndPositionPtrContainer::iterator worstHit;
0458     for (auto it = proto_segment.begin(); it != proto_segment.end(); ++it) {
0459       const float dev = getHitSegChi2(fit, *(*it)->rh);
0460       if (dev < maxDev)
0461         continue;
0462       maxDev = dev;
0463       worstHit = it;
0464     }
0465     LogDebug("GE0SegAlgoRU") << "[GE0SegAlgoRU::pruneBadHits] pruning one hit-> layer: " << (*worstHit)->layer
0466                              << " eta: " << (*worstHit)->gp.eta() << " phi: " << (*worstHit)->gp.phi()
0467                              << " old chi2/dof: " << fit->chi2() / fit->ndof() << std::endl;
0468     proto_segment.erase(worstHit);
0469     fit = makeFit(proto_segment);
0470   }
0471 }
0472 
0473 float GE0SegAlgoRU::getHitSegChi2(const std::unique_ptr<MuonSegFit>& fit, const GEMRecHit& hit) const {
0474   const auto lp = hit.localPosition();
0475   const auto le = hit.localPositionError();
0476   const float du = fit->xdev(lp.x(), lp.z());
0477   const float dv = fit->ydev(lp.y(), lp.z());
0478 
0479   ROOT::Math::SMatrix<double, 2, 2, ROOT::Math::MatRepSym<double, 2> > IC;
0480   IC(0, 0) = le.xx();
0481   IC(1, 0) = le.xy();
0482   IC(1, 1) = le.yy();
0483 
0484   // Invert covariance matrix
0485   IC.Invert();
0486   return du * du * IC(0, 0) + 2. * du * dv * IC(0, 1) + dv * dv * IC(1, 1);
0487 }
0488 
0489 bool GE0SegAlgoRU::hasHitOnLayer(const HitAndPositionPtrContainer& proto_segment, const unsigned int layer) const {
0490   for (const auto* h : proto_segment)
0491     if (h->layer == layer)
0492       return true;
0493   return false;
0494 }
0495 
0496 void GE0SegAlgoRU::compareProtoSegment(std::unique_ptr<MuonSegFit>& current_fit,
0497                                        HitAndPositionPtrContainer& current_proto_segment,
0498                                        const HitAndPosition& new_hit) const {
0499   const HitAndPosition* old_hit = nullptr;
0500   HitAndPositionPtrContainer new_proto_segment = current_proto_segment;
0501 
0502   for (auto it = new_proto_segment.begin(); it != new_proto_segment.end();) {
0503     if ((*it)->layer == new_hit.layer) {
0504       old_hit = *it;
0505       it = new_proto_segment.erase(it);
0506     } else {
0507       ++it;
0508     }
0509   }
0510   if (old_hit == nullptr)
0511     return;
0512   auto new_fit = addHit(new_proto_segment, new_hit);
0513 
0514   //If on the same strip but different BX choose the closest
0515   bool useNew = false;
0516   if (old_hit->lp == new_hit.lp) {
0517     float avgbx = 0;
0518     for (const auto* h : current_proto_segment)
0519       if (old_hit != h)
0520         avgbx += h->rh->BunchX();
0521     avgbx /= float(current_proto_segment.size() - 1);
0522     if (std::abs(avgbx - new_hit.rh->BunchX()) < std::abs(avgbx - old_hit->rh->BunchX()))
0523       useNew = true;
0524   }  //otherwise base it on chi2
0525   else if (new_fit->chi2() < current_fit->chi2())
0526     useNew = true;
0527 
0528   if (useNew) {
0529     current_proto_segment = new_proto_segment;
0530     current_fit = std::move(new_fit);
0531   }
0532 }
0533 
0534 void GE0SegAlgoRU::increaseProtoSegment(const float maxChi2,
0535                                         std::unique_ptr<MuonSegFit>& current_fit,
0536                                         HitAndPositionPtrContainer& current_proto_segment,
0537                                         const HitAndPosition& new_hit) const {
0538   HitAndPositionPtrContainer new_proto_segment = current_proto_segment;
0539   auto new_fit = addHit(new_proto_segment, new_hit);
0540   if (new_fit->chi2() / new_fit->ndof() < maxChi2) {
0541     current_proto_segment = new_proto_segment;
0542     current_fit = std::move(new_fit);
0543   }
0544 }
[ 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!