File indexing completed on 2024-09-07 04:37:39
0001
0002
0003 #include "RecoLocalMuon/CSCRecHitD/src/CSCHitFromStripOnly.h"
0004 #include "RecoLocalMuon/CSCRecHitD/src/CSCStripData.h"
0005 #include "RecoLocalMuon/CSCRecHitD/src/CSCStripHitData.h"
0006 #include "RecoLocalMuon/CSCRecHitD/src/CSCStripHit.h"
0007 #include "RecoLocalMuon/CSCRecHitD/src/CSCPedestalChoice.h"
0008
0009 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0010
0011 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0012 #include "Geometry/CSCGeometry/interface/CSCChamberSpecs.h"
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016
0017 #include <algorithm>
0018 #include <string>
0019 #include <vector>
0020
0021
0022 CSCHitFromStripOnly::CSCHitFromStripOnly(const edm::ParameterSet& ps)
0023 : recoConditions_(nullptr), calcped_(nullptr), ganged_(false) {
0024 useCalib = ps.getParameter<bool>("CSCUseCalibrations");
0025 bool useStaticPedestals = ps.getParameter<bool>("CSCUseStaticPedestals");
0026 int noOfTimeBinsForDynamicPed = ps.getParameter<int>("CSCNoOfTimeBinsForDynamicPedestal");
0027
0028 theThresholdForAPeak = ps.getParameter<double>("CSCStripPeakThreshold");
0029 theThresholdForCluster = ps.getParameter<double>("CSCStripClusterChargeCut");
0030
0031 LogTrace("CSCRecHit") << "[CSCHitFromStripOnly] CSCUseStaticPedestals = " << useStaticPedestals;
0032 if (!useStaticPedestals)
0033 LogTrace("CSCRecHit") << "[CSCHitFromStripOnly] CSCNoOfTimeBinsForDynamicPedestal = " << noOfTimeBinsForDynamicPed;
0034
0035 if (useStaticPedestals) {
0036 calcped_ = new CSCStaticPedestal();
0037 } else {
0038 if (noOfTimeBinsForDynamicPed == 1) {
0039 calcped_ = new CSCDynamicPedestal1();
0040 } else {
0041 calcped_ = new CSCDynamicPedestal2();
0042 }
0043 }
0044 }
0045
0046 CSCHitFromStripOnly::~CSCHitFromStripOnly() { delete calcped_; }
0047
0048
0049
0050
0051
0052
0053
0054 std::vector<CSCStripHit> CSCHitFromStripOnly::runStrip(const CSCDetId& id,
0055 const CSCLayer* layer,
0056 const CSCStripDigiCollection::Range& rstripd) {
0057 std::vector<CSCStripHit> hitsInLayer;
0058
0059
0060 id_ = id;
0061 layer_ = layer;
0062 nstrips_ = layer->chamber()->specs()->nStrips();
0063
0064 setGanged(false);
0065 if (id_.ring() == 4 && layer_->chamber()->specs()->gangedStrips())
0066 setGanged(true);
0067
0068 LogTrace("CSCHitFromStripOnly") << "[CSCHitFromStripOnly::runStrip] id= " << id_ << " nstrips= " << nstrips_
0069 << " ganged strips? " << ganged();
0070
0071 tmax_cluster = 5;
0072
0073
0074
0075
0076
0077
0078 if (useCalib) {
0079 recoConditions_->stripWeights(id, nstrips_, gainWeight);
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090 }
0091
0092
0093 fillPulseHeights(rstripd);
0094 findMaxima(id);
0095
0096
0097 for (size_t imax = 0; imax != theMaxima.size(); ++imax) {
0098
0099 clusterSize = theClusterSize;
0100 theStrips.clear();
0101 strips_adc.clear();
0102 strips_adcRaw.clear();
0103
0104
0105
0106
0107 float strippos = makeCluster(theMaxima[imax] + 1);
0108
0109
0110
0111
0112
0113 if (tmax_cluster < 3) {
0114 theClosestMaximum.push_back(99);
0115 continue;
0116 }
0117
0118 int maximum_to_left = 99;
0119 int maximum_to_right = 99;
0120 if (imax < theMaxima.size() - 1) {
0121 maximum_to_right = theMaxima.at(imax + 1) - theMaxima.at(imax);
0122 }
0123 if (imax > 0) {
0124 maximum_to_left = theMaxima.at(imax - 1) - theMaxima.at(imax);
0125 }
0126 if (std::abs(maximum_to_right) < std::abs(maximum_to_left)) {
0127 theClosestMaximum.push_back(maximum_to_right);
0128 } else {
0129 theClosestMaximum.push_back(maximum_to_left);
0130 }
0131
0132
0133
0134 bool deadStripL = isDeadStrip(id, theMaxima.at(imax) - 1, nstrips_);
0135 bool deadStripR = isDeadStrip(id, theMaxima.at(imax) + 1, nstrips_);
0136 short int aDeadStrip = 0;
0137 if (!deadStripL && !deadStripR) {
0138 aDeadStrip = 0;
0139 } else if (deadStripL && deadStripR) {
0140 aDeadStrip = 255;
0141 } else {
0142 if (deadStripL) {
0143 aDeadStrip = theMaxima.at(imax) - 1;
0144 } else {
0145 aDeadStrip = theMaxima.at(imax) + 1;
0146 }
0147 }
0148
0149
0150
0151 std::vector<int> theL1AStrips;
0152 for (int ila = 0; ila < (int)theStrips.size(); ila++) {
0153 bool stripMatchCounter = false;
0154 for (auto itl1 = rstripd.first; itl1 != rstripd.second; ++itl1) {
0155 int stripNproto = (*itl1).getStrip();
0156 if (!ganged()) {
0157 if (theStrips[ila] == stripNproto) {
0158 stripMatchCounter = true;
0159 auto sz = (*itl1).getOverlappedSample().size();
0160 int L1AbitOnPlace = 0;
0161 for (auto iBit = 0UL; iBit < sz; iBit++) {
0162 L1AbitOnPlace |= ((*itl1).getL1APhase(iBit) << (15 - iBit));
0163 }
0164 theL1AStrips.push_back(theStrips[ila] | L1AbitOnPlace);
0165 }
0166 } else {
0167 for (int tripl = 0; tripl < 3; ++tripl) {
0168 if (theStrips[ila] == (stripNproto + tripl * 16)) {
0169 stripMatchCounter = true;
0170 auto sz = (*itl1).getOverlappedSample().size();
0171 int L1AbitOnPlace = 0;
0172 for (auto iBit = 0UL; iBit < sz; iBit++) {
0173 L1AbitOnPlace |= ((*itl1).getL1APhase(iBit) << (15 - iBit));
0174 }
0175 theL1AStrips.push_back(theStrips[ila] | L1AbitOnPlace);
0176 }
0177 }
0178 }
0179 }
0180 if (!stripMatchCounter) {
0181 theL1AStrips.push_back(theStrips[ila]);
0182 }
0183 }
0184
0185
0186 CSCStripHit striphit(id,
0187 strippos,
0188 tmax_cluster,
0189 theL1AStrips,
0190 strips_adc,
0191 strips_adcRaw,
0192 theConsecutiveStrips.at(imax),
0193 theClosestMaximum.at(imax),
0194 aDeadStrip);
0195 hitsInLayer.push_back(striphit);
0196 }
0197
0198
0199
0200
0201
0202
0203
0204
0205 return hitsInLayer;
0206 }
0207
0208
0209
0210
0211 float CSCHitFromStripOnly::makeCluster(int centerStrip) {
0212 float strippos = -1.;
0213 clusterSize = theClusterSize;
0214 std::vector<CSCStripHitData> stripDataV;
0215
0216
0217
0218
0219 for (int i = 1; i < theClusterSize / 2 + 1; ++i) {
0220 if (centerStrip - i < 1 || centerStrip + i > int(nstrips_)) {
0221
0222 clusterSize = 2 * i - 1;
0223 }
0224 }
0225 for (int i = -clusterSize / 2; i <= clusterSize / 2; ++i) {
0226 CSCStripHitData data = makeStripData(centerStrip, i);
0227 stripDataV.push_back(data);
0228 theStrips.push_back(centerStrip + i);
0229 }
0230 strippos = findHitOnStripPosition(stripDataV, centerStrip);
0231
0232 LogTrace("CSCHitFromStripOnly") << "[CSCHitFromStripOnly::makeCluster] centerStrip= " << centerStrip
0233 << " strippos=" << strippos;
0234
0235 return strippos;
0236 }
0237
0238
0239
0240
0241 CSCStripHitData CSCHitFromStripOnly::makeStripData(int centerStrip, int offset) {
0242 CSCStripHitData prelimData;
0243 int thisStrip = centerStrip + offset;
0244
0245 int tmax = thePulseHeightMap[centerStrip - 1].tmax();
0246 tmax_cluster = tmax;
0247
0248 std::vector<float> adc(4);
0249 std::vector<float> adcRaw(4);
0250
0251
0252
0253 int istart = tmax - 1;
0254 int istop = std::min(tmax + 2, 7);
0255 adc[3] = 0.1;
0256
0257 if (tmax > 2 && tmax < 7) {
0258 int ibin = thisStrip - 1;
0259 if (thePulseHeightMap[ibin].valid()) {
0260 std::copy(
0261 thePulseHeightMap[ibin].ph().begin() + istart, thePulseHeightMap[ibin].ph().begin() + istop + 1, adc.begin());
0262
0263 std::copy(thePulseHeightMap[ibin].phRaw().begin() + istart,
0264 thePulseHeightMap[ibin].phRaw().begin() + istop + 1,
0265 adcRaw.begin());
0266 }
0267 } else {
0268 adc[0] = 0.1;
0269 adc[1] = 0.1;
0270 adc[2] = 0.1;
0271 adc[3] = 0.1;
0272 adcRaw = adc;
0273 LogTrace("CSCRecHit") << "[CSCHitFromStripOnly::makeStripData] Tmax out of range: contact CSC expert!";
0274 }
0275
0276 if (offset == 0) {
0277 prelimData = CSCStripHitData(thisStrip, tmax_cluster, adcRaw, adc);
0278 } else {
0279 int sign = offset > 0 ? 1 : -1;
0280
0281
0282 for (int i = 1; i <= clusterSize / 2; ++i) {
0283
0284 int testStrip = thisStrip + sign * i;
0285 std::vector<int>::iterator otherMax = find(theMaxima.begin(), theMaxima.end(), testStrip - 1);
0286
0287
0288 if (otherMax == theMaxima.end()) {
0289 prelimData = CSCStripHitData(thisStrip, tmax_cluster, adcRaw, adc);
0290 } else {
0291
0292 std::vector<float> adc1(4);
0293 std::vector<float> adcRaw1(4);
0294 std::vector<float> adc2(4);
0295 std::vector<float> adcRaw2(4);
0296
0297 adc1[3] = 0.1;
0298 adc2[3] = 0.1;
0299 adcRaw1[3] = 0.1;
0300 adcRaw2[3] = 0.1;
0301
0302
0303 if (tmax > 2 && tmax < 7) {
0304 int ibin = testStrip - 1;
0305 int jbin = centerStrip - 1;
0306 if (thePulseHeightMap[ibin].valid()) {
0307 std::copy(thePulseHeightMap[ibin].ph().begin() + istart,
0308 thePulseHeightMap[ibin].ph().begin() + istop + 1,
0309 adc1.begin());
0310 std::copy(thePulseHeightMap[ibin].phRaw().begin() + istart,
0311 thePulseHeightMap[ibin].phRaw().begin() + istop + 1,
0312 adcRaw1.begin());
0313 }
0314
0315 if (thePulseHeightMap[jbin].valid()) {
0316 std::copy(thePulseHeightMap[jbin].ph().begin() + istart,
0317 thePulseHeightMap[jbin].ph().begin() + istop + 1,
0318 adc2.begin());
0319
0320 std::copy(thePulseHeightMap[jbin].phRaw().begin() + istart,
0321 thePulseHeightMap[jbin].phRaw().begin() + istop + 1,
0322 adcRaw2.begin());
0323 }
0324 } else {
0325 adc1.assign(4, 0.1);
0326 adcRaw1 = adc1;
0327 adc2.assign(4, 0.1);
0328 adcRaw2 = adc2;
0329 }
0330
0331
0332
0333
0334 for (size_t k = 0; k < 4; ++k) {
0335 if (adc1[k] > 0 && adc2[k] > 0)
0336 adc[k] = adc[k] * adc2[k] / (adc1[k] + adc2[k]);
0337 if (adcRaw1[k] > 0 && adcRaw2[k] > 0)
0338 adcRaw[k] = adcRaw[k] * adcRaw2[k] / (adcRaw1[k] + adcRaw2[k]);
0339 }
0340 prelimData = CSCStripHitData(thisStrip, tmax_cluster, adcRaw, adc);
0341 }
0342 }
0343 }
0344 return prelimData;
0345 }
0346
0347
0348
0349
0350 void CSCHitFromStripOnly::fillPulseHeights(const CSCStripDigiCollection::Range& rstripd) {
0351
0352
0353
0354 for (auto& ph : thePulseHeightMap)
0355 ph.reset();
0356
0357
0358 for (CSCStripDigiCollection::const_iterator it = rstripd.first; it != rstripd.second; ++it) {
0359 int thisChannel = (*it).getStrip();
0360 auto& stripData = thePulseHeightMap[thisChannel - 1];
0361 auto& scaRaw = stripData.phRaw_;
0362 auto& sca = stripData.ph_;
0363
0364 auto const& scaOri = (*it).getADCCounts();
0365 assert(scaOri.size() == 8);
0366
0367 std::copy(scaOri.begin(), scaOri.end(), scaRaw.begin());
0368 std::copy(scaRaw.begin(), scaRaw.end(), sca.begin());
0369
0370
0371 int tmax = std::max_element(sca.begin(), sca.end()) - sca.begin();
0372
0373
0374 float ped = calcped_->pedestal(sca, recoConditions_, id_, thisChannel);
0375
0376
0377 std::for_each(sca.begin(), sca.end(), CSCSubtractPedestal(ped));
0378 std::for_each(scaRaw.begin(), scaRaw.end(), CSCSubtractPedestal(ped));
0379
0380
0381 float phmax = 0.f;
0382 if (tmax > 2 && tmax < 7) {
0383 phmax = sca[tmax];
0384 }
0385 stripData.phmax_ = phmax;
0386 stripData.tmax_ = tmax;
0387
0388
0389
0390
0391
0392
0393 if (useCalib)
0394 stripData *= gainWeight[thisChannel - 1];
0395
0396
0397 if (ganged()) {
0398 for (int j = 1; j < 3; ++j) {
0399 thePulseHeightMap[thisChannel - 1 + 16 * j] = stripData;
0400 }
0401 }
0402 }
0403 }
0404
0405
0406
0407
0408
0409
0410
0411 void CSCHitFromStripOnly::findMaxima(const CSCDetId& id) {
0412 theMaxima.clear();
0413 theConsecutiveStrips.clear();
0414 theClosestMaximum.clear();
0415 for (size_t i = 0; i != thePulseHeightMap.size(); ++i) {
0416
0417 float heightCluster = 0.;
0418
0419 bool maximumFound = false;
0420
0421 if (!isDeadStrip(id, i + 1, nstrips_)) {
0422 if (i == 0) {
0423 heightCluster = thePulseHeightMap[i].phmax() + thePulseHeightMap[i + 1].phmax();
0424
0425 if (thePulseHeightMap[i].phmax() >= thePulseHeightMap[i + 1].phmax() && isPeakOK(i, heightCluster)) {
0426 maximumFound = true;
0427 }
0428
0429 } else if (i == thePulseHeightMap.size() - 1) {
0430 heightCluster = thePulseHeightMap[i - 1].phmax() + thePulseHeightMap[i].phmax();
0431
0432 if (thePulseHeightMap[i].phmax() > thePulseHeightMap[i - 1].phmax() && isPeakOK(i, heightCluster)) {
0433 maximumFound = true;
0434 }
0435
0436 } else {
0437 heightCluster =
0438 thePulseHeightMap[i - 1].phmax() + thePulseHeightMap[i].phmax() + thePulseHeightMap[i + 1].phmax();
0439
0440 if (thePulseHeightMap[i].phmax() > thePulseHeightMap[i - 1].phmax() &&
0441 thePulseHeightMap[i].phmax() >= thePulseHeightMap[i + 1].phmax() && isPeakOK(i, heightCluster)) {
0442 maximumFound = true;
0443 }
0444 }
0445 }
0446
0447 if (maximumFound) {
0448 int numberOfConsecutiveStrips = 1;
0449 float testThreshold = 10.;
0450
0451 int j = 0;
0452 for (int l = 0; l < 8; ++l) {
0453 if (j < 0)
0454 edm::LogWarning("FailedStripCountingWrongConsecutiveStripNumber")
0455 << "This should never occur!!! Contact CSC expert!";
0456 ++j;
0457 bool signalPresent = false;
0458 for (int k = 0; k < 2; ++k) {
0459 j *= -1;
0460 int anotherConsecutiveStrip = i + j;
0461 if (anotherConsecutiveStrip >= 0 && anotherConsecutiveStrip < int(thePulseHeightMap.size())) {
0462 if (thePulseHeightMap[anotherConsecutiveStrip].phmax() > testThreshold) {
0463 ++numberOfConsecutiveStrips;
0464 signalPresent = true;
0465 }
0466 }
0467 }
0468 if (!signalPresent) {
0469 break;
0470 }
0471 }
0472
0473 bool additional_maxima_found = false;
0474
0475
0476
0477
0478
0479 if (i > 2 && i + 3 < thePulseHeightMap.size() && numberOfConsecutiveStrips > 3) {
0480
0481
0482 if (((thePulseHeightMap[i + 1].phmax() >= thePulseHeightMap[i - 1].phmax() &&
0483 thePulseHeightMap[i + 1].phmax() >= thePulseHeightMap[i - 2].phmax() &&
0484 thePulseHeightMap[i + 2].phmax() <= thePulseHeightMap[i - 2].phmax()) ||
0485 (thePulseHeightMap[i + 1].phmax() <= thePulseHeightMap[i - 1].phmax() &&
0486 thePulseHeightMap[i + 1].phmax() <= thePulseHeightMap[i - 2].phmax())) &&
0487
0488 thePulseHeightMap[i - 1].phmax() >= thePulseHeightMap[i - 2].phmax() &&
0489
0490 thePulseHeightMap[i - 2].phmax() > 20) {
0491 additional_maxima_found = true;
0492 theMaxima.push_back(i - 2);
0493
0494 theConsecutiveStrips.push_back(numberOfConsecutiveStrips);
0495 theMaxima.push_back(i);
0496
0497 theConsecutiveStrips.push_back(numberOfConsecutiveStrips);
0498
0499 }
0500
0501
0502 if (((thePulseHeightMap[i + 1].phmax() >= thePulseHeightMap[i - 1].phmax() &&
0503 thePulseHeightMap[i + 2].phmax() >= thePulseHeightMap[i - 1].phmax()) ||
0504 (thePulseHeightMap[i + 1].phmax() <= thePulseHeightMap[i - 1].phmax() &&
0505 thePulseHeightMap[i + 2].phmax() <= thePulseHeightMap[i - 1].phmax() &&
0506 thePulseHeightMap[i + 2].phmax() >= thePulseHeightMap[i - 2].phmax())) &&
0507
0508 thePulseHeightMap[i + 1].phmax() >= thePulseHeightMap[i + 2].phmax() &&
0509
0510 thePulseHeightMap[i + 2].phmax() > 20) {
0511 additional_maxima_found = true;
0512 theMaxima.push_back(i);
0513
0514 theConsecutiveStrips.push_back(numberOfConsecutiveStrips);
0515 theMaxima.push_back(i + 2);
0516
0517 theConsecutiveStrips.push_back(numberOfConsecutiveStrips);
0518
0519 }
0520
0521
0522 if (!additional_maxima_found) {
0523 theMaxima.push_back(i);
0524 theConsecutiveStrips.push_back(numberOfConsecutiveStrips);
0525 }
0526 } else {
0527 theMaxima.push_back(i);
0528 theConsecutiveStrips.push_back(numberOfConsecutiveStrips);
0529 }
0530 }
0531 }
0532 }
0533
0534 bool CSCHitFromStripOnly::isPeakOK(int iStrip, float heightCluster) {
0535 int i = iStrip;
0536 bool peakOK = (thePulseHeightMap[i].phmax() > theThresholdForAPeak && heightCluster > theThresholdForCluster &&
0537
0538
0539 thePulseHeightMap[i].tmax() > 2 && thePulseHeightMap[i].tmax() < 7);
0540 return peakOK;
0541 }
0542
0543
0544
0545
0546 float CSCHitFromStripOnly::findHitOnStripPosition(const std::vector<CSCStripHitData>& data, const int& centerStrip) {
0547 float strippos = -1.;
0548
0549 if (data.empty())
0550 return strippos;
0551
0552
0553
0554
0555 int biggestStrip = max_element(data.begin(), data.end()) - data.begin();
0556 strippos = data[biggestStrip].strip() * 1.;
0557
0558
0559
0560 float sum = 0.;
0561 float sum_w = 0.;
0562
0563
0564
0565
0566 for (size_t i = 0; i != data.size(); ++i) {
0567 auto const& w = data[i].ph();
0568 auto const& wRaw = data[i].phRaw();
0569
0570
0571
0572
0573
0574
0575
0576
0577
0578
0579 std::copy(w.begin(), w.end(), std::back_inserter(strips_adc));
0580 std::copy(wRaw.begin(), wRaw.end(), std::back_inserter(strips_adcRaw));
0581
0582 if (data[i].strip() < 1) {
0583 LogTrace("CSCRecHit") << "[CSCHitFromStripOnly::findHitOnStripPosition] problem in indexing of strip, strip= "
0584 << data[i].strip();
0585 }
0586 sum_w += w[1];
0587 sum += w[1] * data[i].strip();
0588 }
0589
0590 if (sum_w > 0.)
0591 strippos = sum / sum_w;
0592
0593 return strippos;
0594 }
0595
0596 bool CSCHitFromStripOnly::isNearDeadStrip(const CSCDetId& id, int centralStrip, int nstrips) {
0597
0598
0599 return recoConditions_->nearBadStrip(id, centralStrip, nstrips);
0600 }
0601
0602 bool CSCHitFromStripOnly::isDeadStrip(const CSCDetId& id, int centralStrip, int nstrips) {
0603 return recoConditions_->badStrip(id, centralStrip, nstrips);
0604 }
0605
0606
0607 const int CSCHitFromStripOnly::theClusterSize;