Back to home page

Project CMSSW displayed by LXR

 
 

    


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

0001 /**
0002  * \file CSCSegAlgoShowering.cc
0003  *
0004  *  Last update: 17.02.2015
0005  *
0006  */
0007 
0008 #include "RecoLocalMuon/CSCSegment/src/CSCSegAlgoShowering.h"
0009 #include "RecoLocalMuon/CSCSegment/src/CSCSegFit.h"
0010 
0011 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0012 #include <Geometry/CSCGeometry/interface/CSCChamber.h>
0013 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0014 
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 
0018 #include <algorithm>
0019 #include <cmath>
0020 #include <iostream>
0021 #include <string>
0022 
0023 /* Constructor
0024  *
0025  */
0026 CSCSegAlgoShowering::CSCSegAlgoShowering(const edm::ParameterSet& ps) : sfit_(nullptr) {
0027   //  debug                  = ps.getUntrackedParameter<bool>("CSCSegmentDebug");
0028   dRPhiFineMax = ps.getParameter<double>("dRPhiFineMax");
0029   dPhiFineMax = ps.getParameter<double>("dPhiFineMax");
0030   tanThetaMax = ps.getParameter<double>("tanThetaMax");
0031   tanPhiMax = ps.getParameter<double>("tanPhiMax");
0032   maxRatioResidual = ps.getParameter<double>("maxRatioResidualPrune");
0033   //  maxDR                  = ps.getParameter<double>("maxDR");
0034   maxDTheta = ps.getParameter<double>("maxDTheta");
0035   maxDPhi = ps.getParameter<double>("maxDPhi");
0036 }
0037 
0038 /* Destructor:
0039  *
0040  */
0041 CSCSegAlgoShowering::~CSCSegAlgoShowering() {}
0042 
0043 /* showerSeg
0044  *
0045  */
0046 CSCSegment CSCSegAlgoShowering::showerSeg(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
0047   theChamber = aChamber;
0048   // Initialize parameters
0049   std::vector<float> x, y, gz;
0050   std::vector<int> n;
0051 
0052   for (int i = 0; i < 6; ++i) {
0053     x.push_back(0.);
0054     y.push_back(0.);
0055     gz.push_back(0.);
0056     n.push_back(0);
0057   }
0058 
0059   // Loop over hits to find center-of-mass position in each layer
0060   for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); ++it) {
0061     const CSCRecHit2D& hit = (**it);
0062     const CSCDetId id = hit.cscDetId();
0063     int l_id = id.layer();
0064     const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer());
0065     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0066     LocalPoint lp = theChamber->toLocal(gp);
0067 
0068     ++n[l_id - 1];
0069     x[l_id - 1] += lp.x();
0070     y[l_id - 1] += lp.y();
0071     gz[l_id - 1] += gp.z();
0072   }
0073 
0074   // Determine center of mass for each layer and average center of mass for chamber
0075   float avgChamberX = 0.;
0076   float avgChamberY = 0.;
0077   int n_lay = 0;
0078 
0079   for (unsigned i = 0; i < 6; ++i) {
0080     if (n[i] < 1)
0081       continue;
0082 
0083     x[i] = x[i] / n[i];
0084     y[i] = y[i] / n[i];
0085     avgChamberX += x[i];
0086     avgChamberY += y[i];
0087     n_lay++;
0088   }
0089 
0090   if (n_lay > 0) {
0091     avgChamberX = avgChamberX / n_lay;
0092     avgChamberY = avgChamberY / n_lay;
0093   }
0094 
0095   // Create a FakeSegment origin that points back to the IP
0096   // Keep all this in global coordinates until last minute to avoid screwing up +/- signs !
0097 
0098   LocalPoint lpCOM(avgChamberX, avgChamberY, 0.);
0099   GlobalPoint gpCOM = theChamber->toGlobal(lpCOM);
0100   GlobalVector gvCOM(gpCOM.x(), gpCOM.y(), gpCOM.z());
0101 
0102   float Gdxdz = gpCOM.x() / gpCOM.z();
0103   float Gdydz = gpCOM.y() / gpCOM.z();
0104 
0105   // Figure out the intersection of this vector with each layer of the chamber
0106   // by projecting the vector
0107   std::vector<LocalPoint> layerPoints;
0108 
0109   for (size_t i = 0; i != 6; ++i) {
0110     // Get the layer z coordinates in global frame
0111     const CSCLayer* layer = theChamber->layer(i + 1);
0112     LocalPoint temp(0., 0., 0.);
0113     GlobalPoint gp = layer->toGlobal(temp);
0114     float layer_Z = gp.z();
0115 
0116     // Then compute interesection of vector with that plane
0117     float layer_X = Gdxdz * layer_Z;
0118     float layer_Y = Gdydz * layer_Z;
0119     GlobalPoint Gintersect(layer_X, layer_Y, layer_Z);
0120     LocalPoint Lintersect = theChamber->toLocal(Gintersect);
0121 
0122     float layerX = Lintersect.x();
0123     float layerY = Lintersect.y();
0124     float layerZ = Lintersect.z();
0125     LocalPoint layerPoint(layerX, layerY, layerZ);
0126     layerPoints.push_back(layerPoint);
0127   }
0128 
0129   std::vector<float> r_closest;
0130   std::vector<int> id;
0131   for (size_t i = 0; i != 6; ++i) {
0132     id.push_back(-1);
0133     r_closest.push_back(9999.);
0134   }
0135 
0136   int idx = 0;
0137 
0138   // Loop over all hits and find hit closest to com for that layer.
0139   for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); ++it) {
0140     const CSCRecHit2D& hit = (**it);
0141     int layId = hit.cscDetId().layer();
0142 
0143     const CSCLayer* layer = theChamber->layer(layId);
0144     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0145     LocalPoint lp = theChamber->toLocal(gp);
0146 
0147     float d_x = lp.x() - layerPoints[layId - 1].x();
0148     float d_y = lp.y() - layerPoints[layId - 1].y();
0149 
0150     LocalPoint diff(d_x, d_y, 0.);
0151 
0152     if (fabs(diff.mag()) < r_closest[layId - 1]) {
0153       r_closest[layId - 1] = fabs(diff.mag());
0154       id[layId - 1] = idx;
0155     }
0156     ++idx;
0157   }
0158 
0159   // Now fill vector of rechits closest to center of mass:
0160   protoSegment.clear();
0161   idx = 0;
0162 
0163   // Loop over all hits and find hit closest to com for that layer.
0164   for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); ++it) {
0165     const CSCRecHit2D& hit = (**it);
0166     int layId = hit.cscDetId().layer();
0167 
0168     if (idx == id[layId - 1])
0169       protoSegment.push_back(*it);
0170 
0171     ++idx;
0172   }
0173 
0174   // Reorder hits in protosegment
0175   if (gz[0] > 0.) {
0176     if (gz[0] > gz[5]) {
0177       reverse(protoSegment.begin(), protoSegment.end());
0178     }
0179   } else if (gz[0] < 0.) {
0180     if (gz[0] < gz[5]) {
0181       reverse(protoSegment.begin(), protoSegment.end());
0182     }
0183   }
0184 
0185   // Fit the protosegment
0186   updateParameters();
0187 
0188   // If there is one very bad hit on segment, remove it and refit
0189   if (protoSegment.size() > 4)
0190     pruneFromResidual();
0191 
0192   // If any hit on a layer is closer to segment than original, replace it and refit
0193   for (ChamberHitContainer::const_iterator it = rechits.begin(); it != rechits.end(); it++) {
0194     const CSCRecHit2D* h = *it;
0195     int layer = h->cscDetId().layer();
0196     if (isHitNearSegment(h))
0197       compareProtoSegment(h, layer);
0198   }
0199 
0200   // Check again for a bad hit, and remove and refit if necessary
0201   if (sfit_->nhits() > 5)
0202     pruneFromResidual();
0203 
0204   // Does the fitted line point to the IP?
0205   // If it doesn't, the algorithm has probably failed i.e. that's life!
0206 
0207   GlobalVector protoGlobalDir = theChamber->toGlobal(sfit_->localdir());
0208   double protoTheta = protoGlobalDir.theta();
0209   double protoPhi = protoGlobalDir.phi();
0210   double simTheta = gvCOM.theta();
0211   double simPhi = gvCOM.phi();
0212 
0213   float dTheta = fabs(protoTheta - simTheta);
0214   float dPhi = fabs(protoPhi - simPhi);
0215   //  float dR = sqrt(dEta*dEta + dPhi*dPhi);
0216 
0217   // Flag the segment with chi2=-1 of the segment isn't pointing toward origin
0218   // i.e. flag that the algorithm has probably just failed (I presume we expect
0219   // a segment to point to the IP if the algorithm is successful!)
0220 
0221   double theFlag = -1.;
0222   if (dTheta > maxDTheta || dPhi > maxDPhi) {
0223   } else {
0224     theFlag = sfit_->chi2();  // If it points to IP, just pass fit chi2 as usual
0225   }
0226 
0227   // Create an actual CSCSegment - retrieve all info from the fit
0228   CSCSegment temp(sfit_->hits(), sfit_->intercept(), sfit_->localdir(), sfit_->covarianceMatrix(), theFlag);
0229   delete sfit_;
0230   sfit_ = nullptr;
0231 
0232   return temp;
0233 }
0234 
0235 /* isHitNearSegment
0236  *
0237  * Compare rechit with expected position from proto_segment
0238  */
0239 bool CSCSegAlgoShowering::isHitNearSegment(const CSCRecHit2D* hit) const {
0240   const CSCLayer* layer = theChamber->layer(hit->cscDetId().layer());
0241 
0242   // hit phi position in global coordinates
0243   GlobalPoint Hgp = layer->toGlobal(hit->localPosition());
0244   double Hphi = Hgp.phi();
0245   if (Hphi < 0.)
0246     Hphi += 2. * M_PI;
0247   LocalPoint Hlp = theChamber->toLocal(Hgp);
0248   double z = Hlp.z();
0249 
0250   double LocalX = sfit_->xfit(z);
0251   double LocalY = sfit_->yfit(z);
0252   LocalPoint Slp(LocalX, LocalY, z);
0253   GlobalPoint Sgp = theChamber->toGlobal(Slp);
0254   double Sphi = Sgp.phi();
0255   if (Sphi < 0.)
0256     Sphi += 2. * M_PI;
0257   double R = sqrt(Sgp.x() * Sgp.x() + Sgp.y() * Sgp.y());
0258 
0259   double deltaPhi = Sphi - Hphi;
0260   if (deltaPhi > 2. * M_PI)
0261     deltaPhi -= 2. * M_PI;
0262   if (deltaPhi < -2. * M_PI)
0263     deltaPhi += 2. * M_PI;
0264   if (deltaPhi < 0.)
0265     deltaPhi = -deltaPhi;
0266 
0267   double RdeltaPhi = R * deltaPhi;
0268 
0269   if (RdeltaPhi < dRPhiFineMax && deltaPhi < dPhiFineMax)
0270     return true;
0271 
0272   return false;
0273 }
0274 
0275 /* Method addHit
0276  *
0277  * Test if can add hit to proto segment. If so, try to add it.
0278  *
0279  */
0280 bool CSCSegAlgoShowering::addHit(const CSCRecHit2D* aHit, int layer) {
0281   // Return true if hit was added successfully
0282   // Return false if there is already a hit on the same layer
0283 
0284   bool ok = true;
0285 
0286   // Test that we are not trying to add the same hit again
0287   for (ChamberHitContainer::const_iterator it = protoSegment.begin(); it != protoSegment.end(); it++)
0288     if (aHit == (*it))
0289       return false;
0290 
0291   protoSegment.push_back(aHit);
0292 
0293   return ok;
0294 }
0295 
0296 /* Method compareProtoSegment
0297  *      
0298  * For hit coming from the same layer of an existing hit within the proto segment
0299  * test if achieve better chi^2 by using this hit than the other
0300  *
0301  */
0302 void CSCSegAlgoShowering::compareProtoSegment(const CSCRecHit2D* h, int layer) {
0303   // Try adding the hit to existing segment, and removing one from the  same layer
0304   ChamberHitContainer::iterator it;
0305   for (it = protoSegment.begin(); it != protoSegment.end();) {
0306     if ((*it)->cscDetId().layer() == layer) {
0307       it = protoSegment.erase(it);
0308     } else {
0309       ++it;
0310     }
0311   }
0312   bool ok = addHit(h, layer);
0313   if (ok) {
0314     CSCSegFit* newfit = new CSCSegFit(theChamber, protoSegment);
0315     newfit->fit();
0316     if (newfit->chi2() > sfit_->chi2()) {
0317       // new fit is worse: revert to old fit
0318       delete newfit;
0319     } else {
0320       // new fit is better
0321       delete sfit_;
0322       sfit_ = newfit;
0323     }
0324   }
0325 }
0326 
0327 // Look for a hit with a large deviation from fit
0328 
0329 void CSCSegAlgoShowering::pruneFromResidual(void) {
0330   //@@ THIS FUNCTION HAS 3 RETURNS PATHS!
0331 
0332   // Only prune if have at least 5 hits
0333   if (protoSegment.size() < 5)
0334     return;
0335 
0336   // Now Study residuals
0337 
0338   float maxResidual = 0.;
0339   float sumResidual = 0.;
0340   int nHits = 0;
0341 
0342   ChamberHitContainer::iterator ih;
0343   ChamberHitContainer::iterator ibad = protoSegment.end();
0344 
0345   for (ih = protoSegment.begin(); ih != protoSegment.end(); ++ih) {
0346     const CSCRecHit2D& hit = (**ih);
0347     const CSCLayer* layer = theChamber->layer(hit.cscDetId().layer());
0348     GlobalPoint gp = layer->toGlobal(hit.localPosition());
0349     LocalPoint lp = theChamber->toLocal(gp);
0350 
0351     double u = lp.x();
0352     double v = lp.y();
0353     double z = lp.z();
0354 
0355     //    double du = segfit->xdev( u, z );
0356     //    double dv = segfit->ydev( v, z );
0357 
0358     float residual = sfit_->Rdev(u, v, z);  // == sqrt(du*du + dv*dv)
0359 
0360     sumResidual += residual;
0361     ++nHits;
0362     if (residual > maxResidual) {
0363       maxResidual = residual;
0364       ibad = ih;
0365     }
0366   }
0367 
0368   float corrAvgResidual = (sumResidual - maxResidual) / (nHits - 1);
0369 
0370   // Keep all hits
0371   if (maxResidual / corrAvgResidual < maxRatioResidual)
0372     return;
0373 
0374   // Drop worst hit
0375   if (ibad != protoSegment.end())
0376     protoSegment.erase(ibad);
0377 
0378   // Make a new fit
0379   updateParameters();
0380 
0381   return;
0382 }
0383 
0384 void CSCSegAlgoShowering::updateParameters(void) {
0385   // Create fit for the hits in the protosegment & run it
0386   delete sfit_;
0387   sfit_ = new CSCSegFit(theChamber, protoSegment);
0388   sfit_->fit();
0389 }