Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 14:24:41

0001 /**
0002  * \file CSCSegAlgoSK.cc
0003  *
0004  */
0005 
0006 #include "CSCSegAlgoSK.h"
0007 #include "CSCSegFit.h"
0008 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0009 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0010 
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 
0014 #include <algorithm>
0015 #include <cmath>
0016 #include <iostream>
0017 #include <string>
0018 
0019 CSCSegAlgoSK::CSCSegAlgoSK(const edm::ParameterSet& ps)
0020     : CSCSegmentAlgorithm(ps), myName("CSCSegAlgoSK"), sfit_(nullptr) {
0021   debugInfo = ps.getUntrackedParameter<bool>("verboseInfo");
0022 
0023   dRPhiMax = ps.getParameter<double>("dRPhiMax");
0024   dPhiMax = ps.getParameter<double>("dPhiMax");
0025   dRPhiFineMax = ps.getParameter<double>("dRPhiFineMax");
0026   dPhiFineMax = ps.getParameter<double>("dPhiFineMax");
0027   chi2Max = ps.getParameter<double>("chi2Max");
0028   wideSeg = ps.getParameter<double>("wideSeg");
0029   minLayersApart = ps.getParameter<int>("minLayersApart");
0030 
0031   LogDebug("CSC") << myName << " has algorithm cuts set to: \n"
0032                   << "--------------------------------------------------------------------\n"
0033                   << "dRPhiMax     = " << dRPhiMax << '\n'
0034                   << "dPhiMax      = " << dPhiMax << '\n'
0035                   << "dRPhiFineMax = " << dRPhiFineMax << '\n'
0036                   << "dPhiFineMax  = " << dPhiFineMax << '\n'
0037                   << "chi2Max      = " << chi2Max << '\n'
0038                   << "wideSeg      = " << wideSeg << '\n'
0039                   << "minLayersApart = " << minLayersApart << std::endl;
0040 }
0041 
0042 std::vector<CSCSegment> CSCSegAlgoSK::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
0043   theChamber = aChamber;
0044   return buildSegments(rechits);
0045 }
0046 
0047 std::vector<CSCSegment> CSCSegAlgoSK::buildSegments(const ChamberHitContainer& urechits) {
0048   LogDebug("CSC") << "*********************************************";
0049   LogDebug("CSC") << "Start segment building in the new chamber: " << theChamber->specs()->chamberTypeName();
0050   LogDebug("CSC") << "*********************************************";
0051 
0052   ChamberHitContainer rechits = urechits;
0053   LayerIndex layerIndex(rechits.size());
0054 
0055   for (unsigned int i = 0; i < rechits.size(); i++) {
0056     layerIndex[i] = rechits[i]->cscDetId().layer();
0057   }
0058 
0059   double z1 = theChamber->layer(1)->position().z();
0060   double z6 = theChamber->layer(6)->position().z();
0061 
0062   if (z1 > 0.) {
0063     if (z1 > z6) {
0064       reverse(layerIndex.begin(), layerIndex.end());
0065       reverse(rechits.begin(), rechits.end());
0066     }
0067   } else if (z1 < 0.) {
0068     if (z1 < z6) {
0069       reverse(layerIndex.begin(), layerIndex.end());
0070       reverse(rechits.begin(), rechits.end());
0071     }
0072   }
0073 
0074   //  if (debugInfo) dumpHits(rechits);
0075 
0076   if (rechits.size() < 2) {
0077     LogDebug("CSC") << myName << ": " << rechits.size() << "     hit(s) in chamber is not enough to build a segment.\n";
0078     return std::vector<CSCSegment>();
0079   }
0080 
0081   // We have at least 2 hits. We intend to try all possible pairs of hits to start
0082   // segment building. 'All possible' means each hit lies on different layers in the chamber.
0083   // BUT... once a hit has been assigned to a segment, we don't consider
0084   // it again.
0085 
0086   // Choose first hit (as close to IP as possible) h1 and
0087   // second hit (as far from IP as possible) h2
0088   // To do this we iterate over hits in the chamber by layer - pick two layers.
0089   // @@ Require the two layers are at least 3 layers apart. May need tuning?
0090   // Then we iterate over hits within each of these layers and pick h1 and h2 from these.
0091   // If they are 'close enough' we build an empty segment.
0092   // Then try adding hits to this segment.
0093 
0094   // Initialize flags that a given hit has been allocated to a segment
0095   BoolContainer used(rechits.size(), false);
0096 
0097   // Define buffer for segments we build
0098   std::vector<CSCSegment> segments;
0099 
0100   // This is going to point to fits to hits, and its content will be used to create a CSCSegment
0101   sfit_ = nullptr;
0102 
0103   ChamberHitContainerCIt ib = rechits.begin();
0104   ChamberHitContainerCIt ie = rechits.end();
0105 
0106   // Possibly allow 2 passes, second widening scale factor for cuts
0107   windowScale = 1.;  // scale factor for cuts
0108 
0109   int npass = (wideSeg > 1.) ? 2 : 1;
0110 
0111   for (int ipass = 0; ipass < npass; ++ipass) {
0112     for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
0113       bool segok = false;
0114       if (used[i1 - ib])
0115         continue;
0116 
0117       int layer1 = layerIndex[i1 - ib];  //(*i1)->cscDetId().layer();
0118       const CSCRecHit2D* h1 = *i1;
0119 
0120       for (ChamberHitContainerCIt i2 = ie - 1; i2 != i1; --i2) {
0121         if (used[i2 - ib])
0122           continue;
0123 
0124         int layer2 = layerIndex[i2 - ib];  //(*i2)->cscDetId().layer();
0125 
0126         if (abs(layer2 - layer1) < minLayersApart)
0127           break;
0128         const CSCRecHit2D* h2 = *i2;
0129 
0130         if (areHitsCloseInLocalX(h1, h2) && areHitsCloseInGlobalPhi(h1, h2)) {
0131           proto_segment.clear();
0132 
0133           const CSCLayer* l1 = theChamber->layer(layer1);
0134           GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0135           const CSCLayer* l2 = theChamber->layer(layer2);
0136           GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0137           LogDebug("CSC") << "start new segment from hits "
0138                           << "h1: " << gp1 << " - h2: " << gp2 << "\n";
0139 
0140           //@@ TRY ADDING A HIT - AND FIT
0141           if (!addHit(h1, layer1)) {
0142             LogDebug("CSC") << "  failed to add hit h1\n";
0143             continue;
0144           }
0145 
0146           if (!addHit(h2, layer2)) {
0147             LogDebug("CSC") << "  failed to add hit h2\n";
0148             continue;
0149           }
0150 
0151           // Can only add hits if already have a segment
0152           if (sfit_)
0153             tryAddingHitsToSegment(rechits, used, layerIndex, i1, i2);
0154 
0155           // Check no. of hits on segment, and if enough flag them as used
0156           // and store the segment
0157           segok = isSegmentGood(rechits);
0158           if (segok) {
0159             flagHitsAsUsed(rechits, used);
0160             // Copy the proto_segment and its properties inside a CSCSegment.
0161             // Then fill the segment vector..
0162 
0163             if (proto_segment.empty()) {
0164               LogDebug("CSC") << "No segment has been found !!!\n";
0165             } else {
0166               // Create an actual CSCSegment - retrieve all info from the fit
0167               CSCSegment temp(
0168                   sfit_->hits(), sfit_->intercept(), sfit_->localdir(), sfit_->covarianceMatrix(), sfit_->chi2());
0169               delete sfit_;
0170               sfit_ = nullptr;
0171               LogDebug("CSC") << "Found a segment !!!\n";
0172               if (debugInfo)
0173                 dumpSegment(temp);
0174               segments.push_back(temp);
0175             }
0176           }
0177         }  //   h1 & h2 close
0178 
0179         if (segok)
0180           break;
0181       }  //  i2
0182     }    //  i1
0183 
0184     if (segments.size() > 1)
0185       break;  // only change window if no segments found
0186 
0187     // Increase cut windows by factor of wideSeg
0188     windowScale = wideSeg;
0189 
0190   }  //  ipass
0191 
0192   // Give the segments to the CSCChamber
0193   return segments;
0194 }
0195 
0196 void CSCSegAlgoSK::tryAddingHitsToSegment(const ChamberHitContainer& rechits,
0197                                           const BoolContainer& used,
0198                                           const LayerIndex& layerIndex,
0199                                           const ChamberHitContainerCIt i1,
0200                                           const ChamberHitContainerCIt i2) {
0201   // Iterate over the layers with hits in the chamber
0202   // Skip the layers containing the segment endpoints
0203   // Test each hit on the other layers to see if it is near the segment
0204   // If it is, see whether there is already a hit on the segment from the same layer
0205   //    - if so, and there are more than 2 hits on the segment, copy the segment,
0206   //      replace the old hit with the new hit. If the new segment chi2 is better
0207   //      then replace the original segment with the new one (by swap)
0208   //    - if not, copy the segment, add the hit. If the new chi2/dof is still satisfactory
0209   //      then replace the original segment with the new one (by swap)
0210 
0211   ChamberHitContainerCIt ib = rechits.begin();
0212   ChamberHitContainerCIt ie = rechits.end();
0213 
0214   for (ChamberHitContainerCIt i = ib; i != ie; ++i) {
0215     if (i == i1 || i == i2 || used[i - ib])
0216       continue;
0217 
0218     int layer = layerIndex[i - ib];
0219     const CSCRecHit2D* h = *i;
0220     if (isHitNearSegment(h)) {
0221       GlobalPoint gp1 = theChamber->layer(layer)->toGlobal(h->localPosition());
0222       LogDebug("CSC") << "    hit at global " << gp1 << " is near segment\n.";
0223 
0224       // Don't consider alternate hits on layers holding the two starting points
0225       if (hasHitOnLayer(layer)) {
0226         if (proto_segment.size() <= 2) {
0227           LogDebug("CSC") << "    " << proto_segment.size() << " hits on segment...skip hit on same layer.\n";
0228           continue;
0229         }
0230         compareProtoSegment(h, layer);
0231       } else
0232         increaseProtoSegment(h, layer);
0233     }  // h & seg close
0234   }    // i
0235 }
0236 
0237 bool CSCSegAlgoSK::areHitsCloseInLocalX(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const {
0238   float deltaX = (h1->localPosition() - h2->localPosition()).x();
0239   LogDebug("CSC") << "    Hits at local x= " << h1->localPosition().x() << ", " << h2->localPosition().x()
0240                   << " have separation= " << deltaX;
0241   return (fabs(deltaX) < (dRPhiMax * windowScale)) ? true : false;  // +v
0242 }
0243 
0244 bool CSCSegAlgoSK::areHitsCloseInGlobalPhi(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const {
0245   const CSCLayer* l1 = theChamber->layer(h1->cscDetId().layer());
0246   GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0247   const CSCLayer* l2 = theChamber->layer(h2->cscDetId().layer());
0248   GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0249 
0250   float h1p = gp1.phi();
0251   float h2p = gp2.phi();
0252   float dphi12 = h1p - h2p;
0253 
0254   // Into range [-pi, pi) (phi() returns values in this range)
0255   if (dphi12 < -M_PI)
0256     dphi12 += 2. * M_PI;
0257   if (dphi12 > M_PI)
0258     dphi12 -= 2. * M_PI;
0259   LogDebug("CSC") << "    Hits at global phi= " << h1p << ", " << h2p << " have separation= " << dphi12;
0260   return (fabs(dphi12) < (dPhiMax * windowScale)) ? true : false;  // +v
0261 }
0262 
0263 bool CSCSegAlgoSK::isHitNearSegment(const CSCRecHit2D* h) const {
0264   // Is hit near segment?
0265   // Requires deltaphi and rxy*deltaphi within ranges specified
0266   // in parameter set, where rxy=sqrt(x**2+y**2) of hit itself.
0267   // Note that to make intuitive cuts on delta(phi) one must work in
0268   // phi range (-pi, +pi] not [0, 2pi)
0269 
0270   const CSCLayer* l1 = theChamber->layer(h->cscDetId().layer());
0271   GlobalPoint hp = l1->toGlobal(h->localPosition());
0272 
0273   float hphi = hp.phi();  // in (-pi, +pi]
0274   if (hphi < 0.)
0275     hphi += 2. * M_PI;          // into range [0, 2pi)
0276   float sphi = phiAtZ(hp.z());  // in [0, 2*pi)
0277   float phidif = sphi - hphi;
0278   if (phidif < 0.)
0279     phidif += 2. * M_PI;  // into range [0, 2pi)
0280   if (phidif > M_PI)
0281     phidif -= 2. * M_PI;  // into range (-pi, pi]
0282 
0283   float dRPhi = fabs(phidif) * hp.perp();
0284   LogDebug("CSC") << "    is hit at phi_h= " << hphi << " near segment phi_seg= " << sphi << "? is " << dRPhi << "<"
0285                   << dRPhiFineMax * windowScale << " ? "
0286                   << " and is |" << phidif << "|<" << dPhiFineMax * windowScale << " ?";
0287 
0288   return ((dRPhi < dRPhiFineMax * windowScale) && (fabs(phidif) < dPhiFineMax * windowScale)) ? true : false;  // +v
0289 }
0290 
0291 float CSCSegAlgoSK::phiAtZ(float z) const {
0292   if (!sfit_) {
0293     edm::LogVerbatim("CSCSegment") << "[CSCSegAlgoSK::phiAtZ] Segment fit undefined";
0294     return 0.;
0295   }
0296 
0297   // Returns a phi in [ 0, 2*pi )
0298   const CSCLayer* l1 = theChamber->layer((*(proto_segment.begin()))->cscDetId().layer());
0299   GlobalPoint gp = l1->toGlobal(sfit_->intercept());
0300   GlobalVector gv = l1->toGlobal(sfit_->localdir());
0301 
0302   LogTrace("CSCSegment") << "[CSCSegAlgoSK::phiAtZ] Global intercept = " << gp << ", direction = " << gv;
0303 
0304   float x = gp.x() + (gv.x() / gv.z()) * (z - gp.z());
0305   float y = gp.y() + (gv.y() / gv.z()) * (z - gp.z());
0306   float phi = atan2(y, x);
0307   if (phi < 0.f)
0308     phi += 2. * M_PI;
0309 
0310   return phi;
0311 }
0312 
0313 void CSCSegAlgoSK::dumpHits(const ChamberHitContainer& rechits) const {
0314   // Dump positions of RecHit's in each CSCChamber
0315   ChamberHitContainerCIt it;
0316   edm::LogInfo("CSCSegment") << "CSCChamber rechit dump.\n";
0317   for (it = rechits.begin(); it != rechits.end(); it++) {
0318     const CSCLayer* l1 = theChamber->layer((*it)->cscDetId().layer());
0319     GlobalPoint gp1 = l1->toGlobal((*it)->localPosition());
0320 
0321     edm::LogInfo("CSCSegment") << "Global pos.: " << gp1 << ", phi: " << gp1.phi()
0322                                << ". Local position: " << (*it)->localPosition()
0323                                << ", phi: " << (*it)->localPosition().phi() << ". Layer: " << (*it)->cscDetId().layer()
0324                                << "\n";
0325   }
0326 }
0327 
0328 bool CSCSegAlgoSK::isSegmentGood(const ChamberHitContainer& rechitsInChamber) const {
0329   // If the chamber has 20 hits or fewer, require at least 3 hits on segment
0330   // If the chamber has >20 hits require at least 4 hits
0331 
0332   bool ok = false;
0333 
0334   unsigned int iadd = (rechitsInChamber.size() > 20) ? 1 : 0;
0335 
0336   if (windowScale > 1.)
0337     iadd = 1;
0338 
0339   if (proto_segment.size() >= 3 + iadd)
0340     ok = true;
0341 
0342   return ok;
0343 }
0344 
0345 void CSCSegAlgoSK::flagHitsAsUsed(const ChamberHitContainer& rechitsInChamber, BoolContainer& used) const {
0346   // Flag hits on segment as used
0347   ChamberHitContainerCIt ib = rechitsInChamber.begin();
0348   ChamberHitContainerCIt hi, iu;
0349 
0350   for (hi = proto_segment.begin(); hi != proto_segment.end(); ++hi) {
0351     for (iu = ib; iu != rechitsInChamber.end(); ++iu) {
0352       if (*hi == *iu)
0353         used[iu - ib] = true;
0354     }
0355   }
0356 }
0357 
0358 bool CSCSegAlgoSK::addHit(const CSCRecHit2D* aHit, int layer) {
0359   // Return true if hit was added successfully
0360   // (and then parameters are updated).
0361   // Return false if there is already a hit on the same layer, or insert failed.
0362 
0363   ChamberHitContainer::const_iterator it;
0364 
0365   for (it = proto_segment.begin(); it != proto_segment.end(); it++)
0366     if (((*it)->cscDetId().layer() == layer) && (aHit != (*it)))
0367       return false;
0368 
0369   proto_segment.push_back(aHit);
0370 
0371   // make a fit
0372   updateParameters();
0373   return true;
0374 }
0375 
0376 void CSCSegAlgoSK::updateParameters(void) {
0377   // Delete input CSCSegFit, create a new one and make the fit
0378   delete sfit_;
0379   sfit_ = new CSCSegFit(theChamber, proto_segment);
0380   sfit_->fit();
0381 }
0382 
0383 bool CSCSegAlgoSK::hasHitOnLayer(int layer) const {
0384   // Is there is already a hit on this layer?
0385   ChamberHitContainerCIt it;
0386 
0387   for (it = proto_segment.begin(); it != proto_segment.end(); it++)
0388     if ((*it)->cscDetId().layer() == layer)
0389       return true;
0390 
0391   return false;
0392 }
0393 
0394 bool CSCSegAlgoSK::replaceHit(const CSCRecHit2D* h, int layer) {
0395   // replace a hit from a layer
0396   ChamberHitContainer::iterator it;
0397   for (it = proto_segment.begin(); it != proto_segment.end();) {
0398     if ((*it)->cscDetId().layer() == layer)
0399       it = proto_segment.erase(it);
0400     else
0401       ++it;
0402   }
0403 
0404   return addHit(h, layer);
0405 }
0406 
0407 void CSCSegAlgoSK::compareProtoSegment(const CSCRecHit2D* h, int layer) {
0408   // Copy the input CSCSegFit
0409   CSCSegFit* oldfit = new CSCSegFit(*sfit_);
0410 
0411   // May create a new fit
0412   bool ok = replaceHit(h, layer);
0413 
0414   if (ok) {
0415     LogDebug("CSCSegment") << "    hit in same layer as a hit on segment; try replacing old one..."
0416                            << " chi2 new: " << sfit_->chi2() << " old: " << oldfit->chi2() << "\n";
0417   }
0418 
0419   if ((sfit_->chi2() < oldfit->chi2()) && ok) {
0420     LogDebug("CSC") << "    segment with replaced hit is better.\n";
0421     delete oldfit;  // new fit is better
0422   } else {
0423     // keep original fit
0424     delete sfit_;    // now the new fit
0425     sfit_ = oldfit;  // reset to the original input fit
0426   }
0427 }
0428 
0429 void CSCSegAlgoSK::increaseProtoSegment(const CSCRecHit2D* h, int layer) {
0430   // Copy input fit
0431   CSCSegFit* oldfit = new CSCSegFit(*sfit_);
0432 
0433   // Creates a new fit
0434   bool ok = addHit(h, layer);
0435 
0436   if (ok) {
0437     LogDebug("CSCSegment") << "    hit in new layer: added to segment, new chi2: " << sfit_->chi2() << "\n";
0438   }
0439 
0440   //  int ndf = 2*proto_segment.size() - 4;
0441 
0442   //@@ TEST ON ndof<=0 IS JUST TO ACCEPT nhits=2 CASE??
0443   if (ok && ((sfit_->ndof() <= 0) || (sfit_->chi2() / sfit_->ndof() < chi2Max))) {
0444     LogDebug("CSCSegment") << "    segment with added hit is good.\n";
0445     delete oldfit;  // new fit is better
0446   } else {
0447     // reset to original fit
0448     delete sfit_;
0449     sfit_ = oldfit;
0450   }
0451 }
0452 
0453 void CSCSegAlgoSK::dumpSegment(const CSCSegment& seg) const {
0454   edm::LogVerbatim("CSCSegment") << "CSCSegment in " << theChamber->id() << "\nlocal position = " << seg.localPosition()
0455                                  << "\nerror = " << seg.localPositionError()
0456                                  << "\nlocal direction = " << seg.localDirection()
0457                                  << "\nerror =" << seg.localDirectionError() << "\ncovariance matrix"
0458                                  << seg.parametersError() << "chi2/ndf = " << seg.chi2() << "/"
0459                                  << seg.degreesOfFreedom() << "\n#rechits = " << seg.specificRecHits().size()
0460                                  << "\ntime = " << seg.time();
0461 }