File indexing completed on 2021-05-06 01:53:44
0001 #include "FWCore/Utilities/interface/Exception.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0004 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0005 #include "L1Trigger/CSCTriggerPrimitives/plugins/CSCAnodeLCTAnalyzer.h"
0006 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0007
0008 using namespace std;
0009
0010
0011
0012
0013
0014 bool CSCAnodeLCTAnalyzer::debug = true;
0015 bool CSCAnodeLCTAnalyzer::doME1A = true;
0016
0017 vector<CSCAnodeLayerInfo> CSCAnodeLCTAnalyzer::getSimInfo(const CSCALCTDigi& alct,
0018 const CSCDetId& alctId,
0019 const CSCWireDigiCollection* wiredc,
0020 const edm::PSimHitContainer* allSimHits) {
0021
0022
0023
0024
0025
0026 vector<CSCAnodeLayerInfo> alctInfo = lctDigis(alct, alctId, wiredc);
0027
0028
0029 if (alctInfo.size() > CSCConstants::NUM_LAYERS) {
0030 throw cms::Exception("CSCAnodeLCTAnalyzer") << "+++ Number of CSCAnodeLayerInfo objects, " << alctInfo.size()
0031 << ", exceeds max expected, " << CSCConstants::NUM_LAYERS << " +++\n";
0032 }
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042 vector<CSCAnodeLayerInfo>::iterator pali;
0043 for (pali = alctInfo.begin(); pali != alctInfo.end(); pali++) {
0044 digiSimHitAssociator(*pali, allSimHits);
0045 }
0046
0047 return alctInfo;
0048 }
0049
0050 vector<CSCAnodeLayerInfo> CSCAnodeLCTAnalyzer::lctDigis(const CSCALCTDigi& alct,
0051 const CSCDetId& alctId,
0052 const CSCWireDigiCollection* wiredc) {
0053
0054
0055
0056
0057 CSCAnodeLayerInfo tempInfo;
0058 vector<CSCAnodeLayerInfo> vectInfo;
0059
0060
0061 int alct_pattern = 0;
0062 if (!alct.getAccelerator())
0063 alct_pattern = alct.getCollisionB() + 1;
0064 int alct_keywire = alct.getKeyWG();
0065 int alct_bx = alct.getBX();
0066
0067
0068 int MESelection = (alctId.station() < 3) ? 0 : 1;
0069
0070 if (debug) {
0071 LogDebug("lctDigis") << "\nlctDigis: ALCT keywire = " << alct_keywire << "; alctId:"
0072 << " endcap " << alctId.endcap() << ", station " << alctId.station() << ", ring "
0073 << alctId.ring() << ", chamber " << alctId.chamber();
0074 }
0075
0076 for (int i_layer = 0; i_layer < CSCConstants::NUM_LAYERS; i_layer++) {
0077 map<int, CSCWireDigi> digiMap;
0078
0079 tempInfo.clear();
0080
0081
0082
0083
0084 CSCDetId layerId(alctId.endcap(), alctId.station(), alctId.ring(), alctId.chamber(), i_layer + 1);
0085
0086 preselectDigis(alct_bx, layerId, wiredc, digiMap);
0087
0088
0089
0090
0091 if (doME1A) {
0092 if (alctId.station() == 1 && alctId.ring() == 1) {
0093 CSCDetId layerId_me1a(alctId.endcap(), alctId.station(), 4, alctId.chamber(), i_layer + 1);
0094 preselectDigis(alct_bx, layerId_me1a, wiredc, digiMap);
0095 }
0096 }
0097
0098
0099 for (int i_layer = 0; i_layer < CSCConstants::NUM_LAYERS; ++i_layer) {
0100 for (int i_wire = 0; i_wire < CSCConstants::ALCT_PATTERN_WIDTH; ++i_wire) {
0101
0102 if (CSCPatternBank::alct_pattern_legacy_[alct_pattern][i_layer][i_wire]) {
0103 int wire = alct_keywire + CSCPatternBank::alct_keywire_offset_[MESelection][i_wire];
0104 if (wire >= 0 && wire < CSCConstants::MAX_NUM_WIREGROUPS) {
0105
0106 if (digiMap.count(wire) > 0) {
0107 tempInfo.setId(layerId);
0108 tempInfo.addComponent(digiMap[wire]);
0109 if (debug)
0110 LogDebug("lctDigis") << " Digi on ALCT: wire group " << digiMap[wire].getWireGroup();
0111 }
0112 }
0113 }
0114 }
0115 }
0116
0117
0118 if (tempInfo.getId().layer() != 0) {
0119 vectInfo.push_back(tempInfo);
0120 }
0121 }
0122 return vectInfo;
0123 }
0124
0125 void CSCAnodeLCTAnalyzer::preselectDigis(const int alct_bx,
0126 const CSCDetId& layerId,
0127 const CSCWireDigiCollection* wiredc,
0128 map<int, CSCWireDigi>& digiMap) {
0129
0130
0131
0132
0133 const int fifo_tbins = 16;
0134 const int drift_delay = 2;
0135 const int hit_persist = 6;
0136
0137 const CSCWireDigiCollection::Range rwired = wiredc->get(layerId);
0138 for (CSCWireDigiCollection::const_iterator digiIt = rwired.first; digiIt != rwired.second; ++digiIt) {
0139 if (debug)
0140 LogDebug("lctDigis") << "Wire digi: layer " << layerId.layer() - 1 << (*digiIt);
0141 int bx_time = (*digiIt).getTimeBin();
0142 if (bx_time >= 0 && bx_time < fifo_tbins) {
0143
0144 int latch_bx = alct_bx + drift_delay;
0145 if (bx_time <= latch_bx - hit_persist || bx_time > latch_bx) {
0146 if (debug)
0147 LogDebug("lctDigis") << "Late wire digi: layer " << layerId.layer() - 1 << " " << (*digiIt) << " skipping...";
0148 continue;
0149 }
0150
0151 int i_wire = (*digiIt).getWireGroup() - 1;
0152
0153
0154
0155 if (digiMap.count(i_wire) > 0) {
0156 if (digiMap[i_wire].getTimeBin() > bx_time) {
0157 if (debug) {
0158 LogDebug("lctDigis") << " Replacing good wire digi on wire " << i_wire;
0159 }
0160 digiMap.erase(i_wire);
0161 }
0162 }
0163
0164 digiMap[i_wire] = *digiIt;
0165 if (debug) {
0166 LogDebug("lctDigis") << " Good wire digi: wire group " << i_wire;
0167 }
0168 }
0169 }
0170 }
0171
0172 void CSCAnodeLCTAnalyzer::digiSimHitAssociator(CSCAnodeLayerInfo& info, const edm::PSimHitContainer* allSimHits) {
0173
0174
0175
0176
0177 vector<PSimHit> simHits;
0178
0179 vector<CSCWireDigi> thisLayerDigis = info.getRecDigis();
0180 if (!thisLayerDigis.empty()) {
0181 CSCDetId layerId = info.getId();
0182 bool me11 = (layerId.station() == 1) && (layerId.ring() == 1);
0183
0184
0185 for (edm::PSimHitContainer::const_iterator simHitIt = allSimHits->begin(); simHitIt != allSimHits->end();
0186 simHitIt++) {
0187
0188 CSCDetId hitId = (CSCDetId)(*simHitIt).detUnitId();
0189 if (hitId == layerId)
0190 simHits.push_back(*simHitIt);
0191 if (me11) {
0192 CSCDetId layerId_me1a(layerId.endcap(), layerId.station(), 4, layerId.chamber(), layerId.layer());
0193 if (hitId == layerId_me1a)
0194 simHits.push_back(*simHitIt);
0195 }
0196 }
0197
0198 if (!simHits.empty()) {
0199 ostringstream strstrm;
0200 if (debug) {
0201 strstrm << "\nLayer " << layerId.layer() << " has " << simHits.size() << " SimHit(s); eta value(s) = ";
0202 }
0203
0204
0205 for (vector<CSCWireDigi>::const_iterator prd = thisLayerDigis.begin(); prd != thisLayerDigis.end(); prd++) {
0206 double deltaEtaMin = 999.;
0207 double bestHitEta = 999.;
0208 PSimHit* bestHit = nullptr;
0209
0210 int wiregroup = prd->getWireGroup();
0211 double digiEta = getWGEta(layerId, wiregroup - 1);
0212
0213 const CSCLayer* csclayer = geom_->layer(layerId);
0214 for (vector<PSimHit>::iterator psh = simHits.begin(); psh != simHits.end(); psh++) {
0215
0216 LocalPoint hitLP = psh->localPosition();
0217 GlobalPoint hitGP = csclayer->toGlobal(hitLP);
0218 double hitEta = hitGP.eta();
0219 if (debug)
0220 strstrm << hitEta << " ";
0221
0222 double deltaEta = fabs(hitEta - digiEta);
0223 if (deltaEta < deltaEtaMin) {
0224 deltaEtaMin = deltaEta;
0225 bestHit = &(*psh);
0226 bestHitEta = hitEta;
0227 }
0228 }
0229 if (debug) {
0230 strstrm << "\nDigi eta: " << digiEta << ", closest SimHit eta: " << bestHitEta
0231 << ", particle type: " << bestHit->particleType();
0232
0233 }
0234 info.addComponent(*bestHit);
0235 }
0236 if (debug) {
0237 LogDebug("digiSimHitAssociator") << strstrm.str();
0238 }
0239 }
0240 }
0241 }
0242
0243 int CSCAnodeLCTAnalyzer::nearestWG(const vector<CSCAnodeLayerInfo>& allLayerInfo,
0244 double& closestPhi,
0245 double& closestEta) {
0246
0247
0248
0249
0250
0251 int nearestWG = -999;
0252 PSimHit matchedHit;
0253 bool hit_found = false;
0254 CSCDetId layerId;
0255
0256 vector<CSCAnodeLayerInfo>::const_iterator pli;
0257 for (pli = allLayerInfo.begin(); pli != allLayerInfo.end(); pli++) {
0258
0259 if (pli->getId().layer() == CSCConstants::KEY_ALCT_LAYER) {
0260 vector<PSimHit> thisLayerHits = pli->getSimHits();
0261 if (!thisLayerHits.empty()) {
0262
0263
0264 if (thisLayerHits.size() != 1) {
0265 edm::LogWarning("L1CSCTPEmulatorWrongValues")
0266 << "+++ Warning: " << thisLayerHits.size() << " SimHits in key layer " << CSCConstants::KEY_ALCT_LAYER
0267 << "! +++ \n";
0268 for (unsigned i = 0; i < thisLayerHits.size(); i++) {
0269 edm::LogWarning("L1CSCTPEmulatorWrongValues") << " SimHit # " << i << thisLayerHits[i] << "\n";
0270 }
0271 }
0272 matchedHit = thisLayerHits[0];
0273 layerId = pli->getId();
0274 hit_found = true;
0275 break;
0276 }
0277 }
0278 }
0279
0280 if (!hit_found) {
0281 for (pli = allLayerInfo.begin(); pli != allLayerInfo.end(); pli++) {
0282
0283 if (!(pli->getRecDigis()).empty() && !(pli->getSimHits()).empty()) {
0284
0285 vector<PSimHit> thisLayerHits = pli->getSimHits();
0286 matchedHit = thisLayerHits[0];
0287 layerId = pli->getId();
0288 hit_found = true;
0289 break;
0290 }
0291 }
0292 }
0293
0294
0295 if (hit_found) {
0296 const CSCLayer* csclayer = geom_->layer(layerId);
0297 const CSCLayerGeometry* layerGeom = csclayer->geometry();
0298 int nearestW = layerGeom->nearestWire(matchedHit.localPosition());
0299 nearestWG = layerGeom->wireGroup(nearestW);
0300
0301
0302 nearestWG -= 1;
0303 if (nearestWG < 0 || nearestWG >= CSCConstants::MAX_NUM_WIREGROUPS) {
0304 edm::LogWarning("L1CSCTPEmulatorWrongInput")
0305 << "+++ Warning: nearest wire group, " << nearestWG << ", is not in [0-" << CSCConstants::MAX_NUM_WIREGROUPS
0306 << ") interval +++\n";
0307 }
0308
0309 GlobalPoint thisPoint = csclayer->toGlobal(matchedHit.localPosition());
0310 closestPhi = thisPoint.phi();
0311 closestEta = thisPoint.eta();
0312 ostringstream strstrm;
0313 if (debug)
0314 strstrm << "Matched anode phi: " << closestPhi;
0315 if (closestPhi < 0.) {
0316 closestPhi += 2. * M_PI;
0317 if (debug)
0318 strstrm << " (" << closestPhi << ")";
0319 }
0320 if (debug) {
0321 strstrm << " eta: " << closestEta << " on a layer " << layerId.layer() << " (1-6);"
0322 << " nearest wire group: " << nearestWG;
0323 LogDebug("nearestWG") << strstrm.str();
0324 }
0325 }
0326
0327 return nearestWG;
0328 }
0329
0330 void CSCAnodeLCTAnalyzer::setGeometry(const CSCGeometry* geom) { geom_ = geom; }
0331
0332 double CSCAnodeLCTAnalyzer::getWGEta(const CSCDetId& layerId, const int wiregroup) {
0333
0334 if (wiregroup < 0 || wiregroup >= CSCConstants::MAX_NUM_WIREGROUPS) {
0335 edm::LogWarning("L1CSCTPEmulatorWrongInput") << "+++ Warning: wire group, " << wiregroup << ", is not in [0-"
0336 << CSCConstants::MAX_NUM_WIREGROUPS << ") interval +++\n";
0337 }
0338
0339 const CSCLayer* csclayer = geom_->layer(layerId);
0340 const CSCLayerGeometry* layerGeom = csclayer->geometry();
0341 LocalPoint digiLP = layerGeom->localCenterOfWireGroup(wiregroup + 1);
0342
0343
0344
0345
0346 GlobalPoint digiGP = csclayer->toGlobal(digiLP);
0347 double eta = digiGP.eta();
0348
0349 return eta;
0350 }