File indexing completed on 2024-04-06 12:26:00
0001
0002
0003
0004
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
0024
0025
0026 CSCSegAlgoShowering::CSCSegAlgoShowering(const edm::ParameterSet& ps) : sfit_(nullptr) {
0027
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
0034 maxDTheta = ps.getParameter<double>("maxDTheta");
0035 maxDPhi = ps.getParameter<double>("maxDPhi");
0036 }
0037
0038
0039
0040
0041 CSCSegAlgoShowering::~CSCSegAlgoShowering() {}
0042
0043
0044
0045
0046 CSCSegment CSCSegAlgoShowering::showerSeg(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
0047 theChamber = aChamber;
0048
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
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
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
0096
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
0106
0107 std::vector<LocalPoint> layerPoints;
0108
0109 for (size_t i = 0; i != 6; ++i) {
0110
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
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
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
0160 protoSegment.clear();
0161 idx = 0;
0162
0163
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
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
0186 updateParameters();
0187
0188
0189 if (protoSegment.size() > 4)
0190 pruneFromResidual();
0191
0192
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
0201 if (sfit_->nhits() > 5)
0202 pruneFromResidual();
0203
0204
0205
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
0216
0217
0218
0219
0220
0221 double theFlag = -1.;
0222 if (dTheta > maxDTheta || dPhi > maxDPhi) {
0223 } else {
0224 theFlag = sfit_->chi2();
0225 }
0226
0227
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
0236
0237
0238
0239 bool CSCSegAlgoShowering::isHitNearSegment(const CSCRecHit2D* hit) const {
0240 const CSCLayer* layer = theChamber->layer(hit->cscDetId().layer());
0241
0242
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
0276
0277
0278
0279
0280 bool CSCSegAlgoShowering::addHit(const CSCRecHit2D* aHit, int layer) {
0281
0282
0283
0284 bool ok = true;
0285
0286
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
0297
0298
0299
0300
0301
0302 void CSCSegAlgoShowering::compareProtoSegment(const CSCRecHit2D* h, int layer) {
0303
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
0318 delete newfit;
0319 } else {
0320
0321 delete sfit_;
0322 sfit_ = newfit;
0323 }
0324 }
0325 }
0326
0327
0328
0329 void CSCSegAlgoShowering::pruneFromResidual(void) {
0330
0331
0332
0333 if (protoSegment.size() < 5)
0334 return;
0335
0336
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
0356
0357
0358 float residual = sfit_->Rdev(u, v, z);
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
0371 if (maxResidual / corrAvgResidual < maxRatioResidual)
0372 return;
0373
0374
0375 if (ibad != protoSegment.end())
0376 protoSegment.erase(ibad);
0377
0378
0379 updateParameters();
0380
0381 return;
0382 }
0383
0384 void CSCSegAlgoShowering::updateParameters(void) {
0385
0386 delete sfit_;
0387 sfit_ = new CSCSegFit(theChamber, protoSegment);
0388 sfit_->fit();
0389 }