File indexing completed on 2024-09-07 04:37:39
0001
0002
0003
0004
0005
0006
0007
0008 #include "CSCSegAlgoTC.h"
0009 #include "CSCSegFit.h"
0010
0011 #include "DataFormats/CSCRecHit/interface/CSCSegment.h"
0012 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0013
0014 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0015
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "CommonTools/Statistics/interface/ChiSquaredProbability.h"
0020
0021 #include <algorithm>
0022 #include <cmath>
0023 #include <iostream>
0024 #include <string>
0025
0026 CSCSegAlgoTC::CSCSegAlgoTC(const edm::ParameterSet& ps)
0027 : CSCSegmentAlgorithm(ps), sfit_(nullptr), myName("CSCSegAlgoTC") {
0028 debugInfo = ps.getUntrackedParameter<bool>("verboseInfo");
0029
0030 dRPhiMax = ps.getParameter<double>("dRPhiMax");
0031 dPhiMax = ps.getParameter<double>("dPhiMax");
0032 dRPhiFineMax = ps.getParameter<double>("dRPhiFineMax");
0033 dPhiFineMax = ps.getParameter<double>("dPhiFineMax");
0034 chi2Max = ps.getParameter<double>("chi2Max");
0035 chi2ndfProbMin = ps.getParameter<double>("chi2ndfProbMin");
0036 minLayersApart = ps.getParameter<int>("minLayersApart");
0037 SegmentSorting = ps.getParameter<int>("SegmentSorting");
0038
0039 LogDebug("CSC") << myName << " has algorithm cuts set to: \n"
0040 << "--------------------------------------------------------------------\n"
0041 << "dRPhiMax = " << dRPhiMax << '\n'
0042 << "dPhiMax = " << dPhiMax << '\n'
0043 << "dRPhiFineMax = " << dRPhiFineMax << '\n'
0044 << "dPhiFineMax = " << dPhiFineMax << '\n'
0045 << "chi2Max = " << chi2Max << '\n'
0046 << "chi2ndfProbMin = " << chi2ndfProbMin << '\n'
0047 << "minLayersApart = " << minLayersApart << '\n'
0048 << "SegmentSorting = " << SegmentSorting << std::endl;
0049 }
0050
0051 std::vector<CSCSegment> CSCSegAlgoTC::run(const CSCChamber* aChamber, const ChamberHitContainer& rechits) {
0052 theChamber = aChamber;
0053
0054 return buildSegments(rechits);
0055 }
0056
0057 std::vector<CSCSegment> CSCSegAlgoTC::buildSegments(const ChamberHitContainer& urechits) {
0058
0059 ChamberHitContainer rechits(urechits);
0060
0061
0062
0063
0064 if (rechits.size() < 2) {
0065 LogDebug("CSC") << myName << ": " << rechits.size() << " hit(s) in chamber is not enough to build a segment.\n";
0066 return std::vector<CSCSegment>();
0067 }
0068
0069 LayerIndex layerIndex(rechits.size());
0070
0071 for (size_t i = 0; i < rechits.size(); ++i) {
0072 short ilay = rechits[i]->cscDetId().layer();
0073 layerIndex[i] = ilay;
0074
0075
0076 }
0077
0078 double z1 = theChamber->layer(1)->position().z();
0079 double z6 = theChamber->layer(6)->position().z();
0080
0081 if (z1 > 0.) {
0082 if (z1 > z6) {
0083 reverse(layerIndex.begin(), layerIndex.end());
0084 reverse(rechits.begin(), rechits.end());
0085 }
0086 } else if (z1 < 0.) {
0087 if (z1 < z6) {
0088 reverse(layerIndex.begin(), layerIndex.end());
0089 reverse(rechits.begin(), rechits.end());
0090 }
0091 }
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115
0116
0117
0118 std::vector<CSCSegment> segments;
0119
0120 sfit_ = nullptr;
0121
0122 ChamberHitContainerCIt ib = rechits.begin();
0123 ChamberHitContainerCIt ie = rechits.end();
0124
0125 for (ChamberHitContainerCIt i1 = ib; i1 != ie; ++i1) {
0126 int layer1 = layerIndex[i1 - ib];
0127
0128 const CSCRecHit2D* h1 = *i1;
0129
0130 for (ChamberHitContainerCIt i2 = ie - 1; i2 != i1; --i2) {
0131 int layer2 = layerIndex[i2 - ib];
0132
0133 if (abs(layer2 - layer1) < minLayersApart)
0134 break;
0135
0136 const CSCRecHit2D* h2 = *i2;
0137
0138 if (areHitsCloseInLocalX(h1, h2) && areHitsCloseInGlobalPhi(h1, h2)) {
0139 proto_segment.clear();
0140
0141 const CSCLayer* l1 = theChamber->layer(h1->cscDetId().layer());
0142 GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0143 const CSCLayer* l2 = theChamber->layer(h2->cscDetId().layer());
0144 GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0145 LogDebug("CSCSegment") << "start new segment from hits "
0146 << "h1: " << gp1 << " - h2: " << gp2 << "\n";
0147
0148
0149
0150 if (!addHit(h1, layer1)) {
0151 LogDebug("CSCSegment") << " failed to add hit h1\n";
0152 continue;
0153 }
0154
0155 if (!addHit(h2, layer2)) {
0156 LogDebug("CSCSegment") << " failed to add hit h2\n";
0157 continue;
0158 }
0159
0160 if (sfit_)
0161 tryAddingHitsToSegment(rechits, i1, i2);
0162
0163
0164 if (proto_segment.empty()) {
0165 LogDebug("CSCSegment") << "No segment found.\n";
0166
0167 } else {
0168
0169
0170 candidates.push_back(sfit_);
0171 LogDebug("CSCSegment") << "Found a segment.\n";
0172
0173 }
0174 }
0175 }
0176 }
0177
0178
0179
0180
0181 pruneTheSegments(rechits);
0182
0183
0184 for (unsigned int i = 0; i < candidates.size(); ++i) {
0185 CSCSegFit* sfit = candidates[i];
0186
0187
0188
0189
0190
0191 CSCSegment temp(sfit->hits(), sfit->intercept(), sfit->localdir(), sfit->covarianceMatrix(), sfit->chi2());
0192 delete sfit;
0193 segments.push_back(temp);
0194 if (debugInfo)
0195 dumpSegment(temp);
0196 }
0197
0198
0199 candidates.clear();
0200 sfit_ = nullptr;
0201
0202
0203
0204 return segments;
0205 }
0206
0207 void CSCSegAlgoTC::tryAddingHitsToSegment(const ChamberHitContainer& rechits,
0208 const ChamberHitContainerCIt i1,
0209 const ChamberHitContainerCIt i2) {
0210
0211
0212
0213
0214
0215
0216
0217
0218
0219
0220 ChamberHitContainerCIt ib = rechits.begin();
0221 ChamberHitContainerCIt ie = rechits.end();
0222
0223 for (ChamberHitContainerCIt i = ib; i != ie; ++i) {
0224 if (i == i1 || i == i2)
0225 continue;
0226
0227 int layer = (*i)->cscDetId().layer();
0228 const CSCRecHit2D* h = *i;
0229
0230 if (isHitNearSegment(h)) {
0231 const CSCLayer* l1 = theChamber->layer(h->cscDetId().layer());
0232 GlobalPoint gp1 = l1->toGlobal(h->localPosition());
0233 LogDebug("CSC") << " hit at global " << gp1 << " is near segment\n.";
0234
0235 if (hasHitOnLayer(layer)) {
0236 if (proto_segment.size() <= 2) {
0237 LogDebug("CSC") << " " << proto_segment.size() << " hits on segment...skip hit on same layer.\n";
0238 continue;
0239 }
0240
0241 compareProtoSegment(h, layer);
0242 } else
0243 increaseProtoSegment(h, layer);
0244 }
0245 }
0246 }
0247
0248 bool CSCSegAlgoTC::addHit(const CSCRecHit2D* aHit, int layer) {
0249
0250
0251
0252
0253 bool ok = true;
0254 ChamberHitContainer::const_iterator it;
0255
0256 for (it = proto_segment.begin(); it != proto_segment.end(); it++)
0257 if (((*it)->cscDetId().layer() == layer) && (aHit != *it))
0258 return false;
0259
0260 if (ok) {
0261 proto_segment.push_back(aHit);
0262 updateParameters();
0263 }
0264 return ok;
0265 }
0266
0267 bool CSCSegAlgoTC::replaceHit(const CSCRecHit2D* h, int layer) {
0268
0269 ChamberHitContainer::iterator it;
0270 for (it = proto_segment.begin(); it != proto_segment.end();) {
0271 if ((*it)->cscDetId().layer() == layer)
0272 it = proto_segment.erase(it);
0273 else
0274 ++it;
0275 }
0276
0277 return addHit(h, layer);
0278 }
0279
0280 void CSCSegAlgoTC::updateParameters(void) {
0281
0282
0283 sfit_ = new CSCSegFit(theChamber, proto_segment);
0284 sfit_->fit();
0285 }
0286
0287 float CSCSegAlgoTC::phiAtZ(float z) const {
0288
0289 const CSCLayer* l1 = theChamber->layer(1);
0290 GlobalPoint gp = l1->toGlobal(sfit_->intercept());
0291 GlobalVector gv = l1->toGlobal(sfit_->localdir());
0292
0293 float x = gp.x() + (gv.x() / gv.z()) * (z - gp.z());
0294 float y = gp.y() + (gv.y() / gv.z()) * (z - gp.z());
0295 float phi = atan2(y, x);
0296 if (phi < 0.f)
0297 phi += 2. * M_PI;
0298
0299 return phi;
0300 }
0301
0302 void CSCSegAlgoTC::compareProtoSegment(const CSCRecHit2D* h, int layer) {
0303
0304 CSCSegFit* oldfit = new CSCSegFit(*sfit_);
0305
0306 bool ok = replaceHit(h, layer);
0307
0308 if (ok) {
0309 LogDebug("CSC") << " hit in same layer as a hit on segment; try replacing old one..."
0310 << " chi2 new: " << sfit_->chi2() << " old: " << oldfit->chi2() << "\n";
0311 }
0312
0313 if (ok && (sfit_->chi2() < oldfit->chi2())) {
0314 LogDebug("CSC") << " segment with replaced hit is better.\n";
0315 delete oldfit;
0316 } else {
0317 delete sfit_;
0318 sfit_ = oldfit;
0319 }
0320 }
0321
0322 void CSCSegAlgoTC::increaseProtoSegment(const CSCRecHit2D* h, int layer) {
0323
0324 CSCSegFit* oldfit = new CSCSegFit(*sfit_);
0325
0326 bool ok = addHit(h, layer);
0327
0328 if (ok) {
0329 LogDebug("CSC") << " hit in new layer: added to segment, new chi2: " << sfit_->chi2() << "\n";
0330 }
0331
0332
0333
0334 if (ok && ((sfit_->chi2() <= 0) || (sfit_->chi2() / sfit_->ndof() < chi2Max))) {
0335 LogDebug("CSC") << " segment with added hit is good.\n";
0336 delete oldfit;
0337 } else {
0338 delete sfit_;
0339 sfit_ = oldfit;
0340 }
0341 }
0342
0343 bool CSCSegAlgoTC::areHitsCloseInLocalX(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const {
0344 float deltaX = (h1->localPosition() - h2->localPosition()).x();
0345 LogDebug("CSC") << " Hits at local x= " << h1->localPosition().x() << ", " << h2->localPosition().x()
0346 << " have separation= " << deltaX;
0347 return (fabs(deltaX) < (dRPhiMax)) ? true : false;
0348 }
0349
0350 bool CSCSegAlgoTC::areHitsCloseInGlobalPhi(const CSCRecHit2D* h1, const CSCRecHit2D* h2) const {
0351 const CSCLayer* l1 = theChamber->layer(h1->cscDetId().layer());
0352 GlobalPoint gp1 = l1->toGlobal(h1->localPosition());
0353 const CSCLayer* l2 = theChamber->layer(h2->cscDetId().layer());
0354 GlobalPoint gp2 = l2->toGlobal(h2->localPosition());
0355
0356 float h1p = gp1.phi();
0357 float h2p = gp2.phi();
0358 float dphi12 = h1p - h2p;
0359
0360
0361 if (dphi12 < -M_PI)
0362 dphi12 += 2. * M_PI;
0363 if (dphi12 > M_PI)
0364 dphi12 -= 2. * M_PI;
0365 LogDebug("CSC") << " Hits at global phi= " << h1p << ", " << h2p << " have separation= " << dphi12;
0366 return (fabs(dphi12) < dPhiMax) ? true : false;
0367 }
0368
0369 bool CSCSegAlgoTC::isHitNearSegment(const CSCRecHit2D* h) const {
0370
0371
0372
0373
0374
0375
0376 const CSCLayer* l1 = theChamber->layer(h->cscDetId().layer());
0377 GlobalPoint hp = l1->toGlobal(h->localPosition());
0378
0379 float hphi = hp.phi();
0380 if (hphi < 0.)
0381 hphi += 2. * M_PI;
0382 float sphi = phiAtZ(hp.z());
0383 float phidif = sphi - hphi;
0384 if (phidif < 0.)
0385 phidif += 2. * M_PI;
0386 if (phidif > M_PI)
0387 phidif -= 2. * M_PI;
0388
0389 float dRPhi = fabs(phidif) * hp.perp();
0390 LogDebug("CSC") << " is hit at phi_h= " << hphi << " near segment phi_seg= " << sphi << "? is " << dRPhi << "<"
0391 << dRPhiFineMax << " ? "
0392 << " and is |" << phidif << "|<" << dPhiFineMax << " ?";
0393
0394 return ((dRPhi < dRPhiFineMax) && (fabs(phidif) < dPhiFineMax)) ? true : false;
0395 }
0396
0397 bool CSCSegAlgoTC::hasHitOnLayer(int layer) const {
0398
0399 ChamberHitContainerCIt it;
0400
0401 for (it = proto_segment.begin(); it != proto_segment.end(); it++)
0402 if ((*it)->cscDetId().layer() == layer)
0403 return true;
0404
0405 return false;
0406 }
0407
0408 void CSCSegAlgoTC::dumpHits(const ChamberHitContainer& rechits) const {
0409
0410 ChamberHitContainerCIt it;
0411
0412 for (it = rechits.begin(); it != rechits.end(); it++) {
0413 const CSCLayer* l1 = theChamber->layer((*it)->cscDetId().layer());
0414 GlobalPoint gp1 = l1->toGlobal((*it)->localPosition());
0415
0416 LogDebug("CSC") << "Global pos.: " << gp1 << ", phi: " << gp1.phi()
0417 << ". Local position: " << (*it)->localPosition() << ", phi: " << (*it)->localPosition().phi()
0418 << ". Layer: " << (*it)->cscDetId().layer() << "\n";
0419 }
0420 }
0421
0422 bool CSCSegAlgoTC::isSegmentGood(std::vector<CSCSegFit*>::iterator seg,
0423 const ChamberHitContainer& rechitsInChamber,
0424 BoolContainer& used) const {
0425
0426
0427
0428
0429
0430
0431
0432
0433 size_t iadd = (rechitsInChamber.size() > 20) ? 1 : 0;
0434
0435 size_t nhits = (*seg)->nhits();
0436
0437 if (nhits < 3 + iadd)
0438 return false;
0439
0440
0441
0442
0443
0444
0445 if (SegmentSorting == 2) {
0446 double chi2t = (*seg)->chi2();
0447 double ndoft = 2 * nhits - 4;
0448 if (chi2t > 0 && ndoft > 0) {
0449 if (ChiSquaredProbability(chi2t, ndoft) < chi2ndfProbMin) {
0450 return false;
0451 }
0452 } else {
0453 return false;
0454 }
0455 }
0456
0457 ChamberHitContainer hits_ = (*seg)->hits();
0458
0459 for (size_t ish = 0; ish < nhits; ++ish) {
0460 ChamberHitContainerCIt ib = rechitsInChamber.begin();
0461
0462 for (ChamberHitContainerCIt ic = ib; ic != rechitsInChamber.end(); ++ic) {
0463 if ((hits_[ish] == (*ic)) && used[ic - ib])
0464 return false;
0465 }
0466 }
0467
0468 return true;
0469 }
0470
0471 void CSCSegAlgoTC::flagHitsAsUsed(std::vector<CSCSegFit*>::iterator seg,
0472 const ChamberHitContainer& rechitsInChamber,
0473 BoolContainer& used) const {
0474
0475
0476 ChamberHitContainerCIt ib = rechitsInChamber.begin();
0477 ChamberHitContainer hits = (*seg)->hits();
0478
0479 for (size_t ish = 0; ish < hits.size(); ++ish) {
0480 for (ChamberHitContainerCIt iu = ib; iu != rechitsInChamber.end(); ++iu)
0481 if (hits[ish] == (*iu))
0482 used[iu - ib] = true;
0483 }
0484 }
0485
0486 void CSCSegAlgoTC::pruneTheSegments(const ChamberHitContainer& rechitsInChamber) {
0487
0488
0489
0490 if (candidates.empty())
0491 return;
0492
0493
0494 BoolContainer used(rechitsInChamber.size(), false);
0495
0496
0497 segmentSort();
0498
0499
0500
0501
0502
0503 std::vector<CSCSegFit*>::iterator is;
0504
0505 for (is = candidates.begin(); is != candidates.end();) {
0506 bool goodSegment = isSegmentGood(is, rechitsInChamber, used);
0507
0508 if (goodSegment) {
0509 LogDebug("CSC") << "Accepting segment: ";
0510 flagHitsAsUsed(is, rechitsInChamber, used);
0511 ++is;
0512 } else {
0513 LogDebug("CSC") << "Rejecting segment: ";
0514 delete *is;
0515 is = candidates.erase(is);
0516 }
0517 }
0518 }
0519
0520 void CSCSegAlgoTC::segmentSort() {
0521
0522
0523 for (size_t i = 0; i < candidates.size() - 1; ++i) {
0524 for (size_t j = i + 1; j < candidates.size(); ++j) {
0525 size_t ni = (candidates[i]->hits()).size();
0526 size_t nj = (candidates[j]->hits()).size();
0527
0528
0529 if (SegmentSorting == 2) {
0530 if (nj > ni) {
0531 CSCSegFit* temp = candidates[j];
0532 candidates[j] = candidates[i];
0533 candidates[i] = temp;
0534 }
0535
0536
0537 if (candidates[i]->chi2() > 0 && candidates[i]->ndof() > 0) {
0538 if (nj == ni && (ChiSquaredProbability(candidates[i]->chi2(), (double)(candidates[i]->ndof())) <
0539 ChiSquaredProbability(candidates[j]->chi2(), (double)(candidates[j]->ndof())))) {
0540 CSCSegFit* temp = candidates[j];
0541 candidates[j] = candidates[i];
0542 candidates[i] = temp;
0543 }
0544 }
0545 } else if (SegmentSorting == 1) {
0546 if ((candidates[i]->chi2() / ni) > (candidates[j]->chi2() / nj)) {
0547 CSCSegFit* temp = candidates[j];
0548 candidates[j] = candidates[i];
0549 candidates[i] = temp;
0550 }
0551 } else {
0552 LogDebug("CSC") << "No valid segment sorting specified. Algorithm misconfigured! \n";
0553 }
0554 }
0555 }
0556 }
0557
0558 void CSCSegAlgoTC::dumpSegment(const CSCSegment& seg) const {
0559 edm::LogVerbatim("CSCSegment") << "CSCSegment in " << theChamber->id() << "\nlocal position = " << seg.localPosition()
0560 << "\nerror = " << seg.localPositionError()
0561 << "\nlocal direction = " << seg.localDirection()
0562 << "\nerror =" << seg.localDirectionError() << "\ncovariance matrix"
0563 << seg.parametersError() << "chi2/ndf = " << seg.chi2() << "/"
0564 << seg.degreesOfFreedom() << "\n#rechits = " << seg.specificRecHits().size()
0565 << "\ntime = " << seg.time();
0566 }