File indexing completed on 2021-02-14 14:14:25
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include "Calibration/Tools/interface/DetIdAssociator.h"
0021
0022 #include <memory>
0023
0024
0025 std::vector<GlobalPoint> HDetIdAssociator::getTrajectory(const FreeTrajectoryState& ftsStart,
0026 const std::vector<GlobalPoint>& surfaces) {
0027 check_setup();
0028 std::vector<GlobalPoint> trajectory;
0029 TrajectoryStateOnSurface tSOSDest;
0030 FreeTrajectoryState ftsCurrent = ftsStart;
0031
0032 for (std::vector<GlobalPoint>::const_iterator surface_iter = surfaces.begin(); surface_iter != surfaces.end();
0033 surface_iter++) {
0034
0035 std::unique_ptr<Cylinder> cylinder =
0036 std::make_unique<Cylinder>(surface_iter->perp(), Surface::PositionType(0, 0, 0), Surface::RotationType());
0037 std::unique_ptr<Plane> forwardEndcap =
0038 std::make_unique<Plane>(Surface::PositionType(0, 0, surface_iter->z()), Surface::RotationType());
0039 std::unique_ptr<Plane> backwardEndcap =
0040 std::make_unique<Plane>(Surface::PositionType(0, 0, -surface_iter->z()), Surface::RotationType());
0041
0042 LogTrace("StartingPoint") << "Propagate from "
0043 << "\n"
0044 << "\tx: " << ftsStart.position().x() << "\n"
0045 << "\ty: " << ftsStart.position().y() << "\n"
0046 << "\tz: " << ftsStart.position().z() << "\n"
0047 << "\tmomentum eta: " << ftsStart.momentum().eta() << "\n"
0048 << "\tmomentum phi: " << ftsStart.momentum().phi() << "\n"
0049 << "\tmomentum: " << ftsStart.momentum().mag() << "\n";
0050
0051 float tanTheta = ftsCurrent.momentum().perp() / ftsCurrent.momentum().z();
0052 float corner = surface_iter->perp() / surface_iter->z();
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066 int ibar = 0;
0067 if (fabs(tanTheta) > corner) {
0068 tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0069
0070 } else if (tanTheta > 0.) {
0071 tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0072 ibar = 1;
0073 } else {
0074 tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0075 ibar = -1;
0076 }
0077
0078
0079
0080 if (!tSOSDest.isValid()) {
0081
0082 if (ibar == 0) {
0083 if (tanTheta < 0)
0084 tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0085 if (tanTheta >= 0)
0086 tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0087 } else {
0088 tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0089 }
0090 } else {
0091
0092 if (abs(ibar) > 0) {
0093 if (tSOSDest.globalPosition().perp() > surface_iter->perp()) {
0094 tSOSDest = ivProp_->propagate(ftsCurrent, *cylinder);
0095 }
0096 } else {
0097 if (tanTheta < 0)
0098 tSOSDest = ivProp_->propagate(ftsCurrent, *forwardEndcap);
0099 if (tanTheta >= 0)
0100 tSOSDest = ivProp_->propagate(ftsCurrent, *backwardEndcap);
0101 }
0102 }
0103
0104
0105
0106
0107
0108
0109
0110
0111
0112
0113
0114
0115 if (!tSOSDest.isValid())
0116 return trajectory;
0117
0118
0119 LogTrace("SuccessfullPropagation") << "Great, I reached something."
0120 << "\n"
0121 << "\tx: " << tSOSDest.freeState()->position().x() << "\n"
0122 << "\ty: " << tSOSDest.freeState()->position().y() << "\n"
0123 << "\tz: " << tSOSDest.freeState()->position().z() << "\n"
0124 << "\teta: " << tSOSDest.freeState()->position().eta() << "\n"
0125 << "\tphi: " << tSOSDest.freeState()->position().phi() << "\n";
0126
0127
0128
0129 GlobalPoint point = tSOSDest.freeState()->position();
0130 point = tSOSDest.freeState()->position();
0131 ftsCurrent = *tSOSDest.freeState();
0132 trajectory.push_back(point);
0133 }
0134 return trajectory;
0135 }
0136
0137
0138 std::set<DetId> HDetIdAssociator::getDetIdsCloseToAPoint(const GlobalPoint& direction, const int idR) {
0139 std::set<DetId> set;
0140 check_setup();
0141 int nDets = 0;
0142 if (!theMap_)
0143 buildMap();
0144 LogTrace("MatchPoint") << "point (eta,phi): " << direction.eta() << "," << direction.phi() << "\n";
0145 int ieta = iEta(direction);
0146 int iphi = iPhi(direction);
0147
0148 LogTrace("MatchPoint") << "(ieta,iphi): " << ieta << "," << iphi << "\n";
0149
0150 if (ieta >= 0 && ieta < nEta_ && iphi >= 0 && iphi < nPhi_) {
0151 set = (*theMap_)[ieta][iphi];
0152 nDets++;
0153 if (idR > 0) {
0154 LogTrace("MatchPoint") << "Add neighbors (ieta,iphi): " << ieta << "," << iphi << "\n";
0155
0156 int maxIEta = ieta + idR;
0157 int minIEta = ieta - idR;
0158 if (maxIEta >= nEta_)
0159 maxIEta = nEta_ - 1;
0160 if (minIEta < 0)
0161 minIEta = 0;
0162 int maxIPhi = iphi + idR;
0163 int minIPhi = iphi - idR;
0164 if (minIPhi < 0) {
0165 minIPhi += nPhi_;
0166 maxIPhi += nPhi_;
0167 }
0168 LogTrace("MatchPoint") << "\tieta (min,max): " << minIEta << "," << maxIEta << "\n";
0169 LogTrace("MatchPoint") << "\tiphi (min,max): " << minIPhi << "," << maxIPhi << "\n";
0170 for (int i = minIEta; i <= maxIEta; i++)
0171 for (int j = minIPhi; j <= maxIPhi; j++) {
0172 if (i == ieta && j == iphi)
0173 continue;
0174 set.insert((*theMap_)[i][j % nPhi_].begin(), (*theMap_)[i][j % nPhi_].end());
0175 nDets++;
0176 }
0177 }
0178 }
0179
0180
0181
0182
0183
0184
0185
0186
0187
0188
0189
0190 return set;
0191 }
0192
0193
0194 int HDetIdAssociator::iEta(const GlobalPoint& point) {
0195
0196 int iEta1 = int(point.eta() / etaBinSize_ + nEta_ / 2);
0197 if (point.eta() > 1.827 && point.eta() <= 1.830)
0198 return iEta1 - 1;
0199 else if (point.eta() > 1.914 && point.eta() <= 1.930)
0200 return iEta1 - 1;
0201 else if (point.eta() > 2.001 && point.eta() <= 2.043)
0202 return iEta1 - 1;
0203 else if (point.eta() > 2.088 && point.eta() <= 2.172)
0204 return iEta1 - 1;
0205 else if (point.eta() > 2.175 && point.eta() <= 2.262)
0206 return iEta1 - 1;
0207 else if (point.eta() > 2.262 && point.eta() <= 2.332)
0208 return iEta1 - 2;
0209 else if (point.eta() > 2.332 && point.eta() <= 2.349)
0210 return iEta1 - 1;
0211 else if (point.eta() > 2.349 && point.eta() <= 2.436)
0212 return iEta1 - 2;
0213 else if (point.eta() > 2.436 && point.eta() <= 2.500)
0214 return iEta1 - 3;
0215 else if (point.eta() > 2.500 && point.eta() <= 2.523)
0216 return iEta1 - 2;
0217 else if (point.eta() > 2.523 && point.eta() <= 2.610)
0218 return iEta1 - 3;
0219 else if (point.eta() > 2.610 && point.eta() <= 2.650)
0220 return iEta1 - 4;
0221 else if (point.eta() > 2.650 && point.eta() <= 2.697)
0222 return iEta1 - 3;
0223 else if (point.eta() > 2.697 && point.eta() <= 2.784)
0224 return iEta1 - 4;
0225 else if (point.eta() > 2.784 && point.eta() <= 2.868)
0226 return iEta1 - 5;
0227 else if (point.eta() > 2.868 && point.eta() <= 2.871)
0228 return iEta1 - 4;
0229 else if (point.eta() > 2.871 && point.eta() <= 2.958)
0230 return iEta1 - 5;
0231 else if (point.eta() > 2.958)
0232 return iEta1 - 6;
0233 else if (point.eta() < -1.827 && point.eta() >= -1.830)
0234 return iEta1 + 1;
0235 else if (point.eta() < -1.914 && point.eta() >= -1.930)
0236 return iEta1 + 1;
0237 else if (point.eta() < -2.001 && point.eta() >= -2.043)
0238 return iEta1 + 1;
0239 else if (point.eta() < -2.088 && point.eta() >= -2.172)
0240 return iEta1 + 1;
0241 else if (point.eta() < -2.175 && point.eta() >= -2.262)
0242 return iEta1 + 1;
0243 else if (point.eta() < -2.262 && point.eta() >= -2.332)
0244 return iEta1 + 2;
0245 else if (point.eta() < -2.332 && point.eta() >= -2.349)
0246 return iEta1 + 1;
0247 else if (point.eta() < -2.349 && point.eta() >= -2.436)
0248 return iEta1 + 2;
0249 else if (point.eta() < -2.436 && point.eta() >= -2.500)
0250 return iEta1 + 3;
0251 else if (point.eta() < -2.500 && point.eta() >= -2.523)
0252 return iEta1 + 2;
0253 else if (point.eta() < -2.523 && point.eta() >= -2.610)
0254 return iEta1 + 3;
0255 else if (point.eta() < -2.610 && point.eta() >= -2.650)
0256 return iEta1 + 4;
0257 else if (point.eta() < -2.650 && point.eta() >= -2.697)
0258 return iEta1 + 3;
0259 else if (point.eta() < -2.697 && point.eta() >= -2.784)
0260 return iEta1 + 4;
0261 else if (point.eta() < -2.784 && point.eta() >= -2.868)
0262 return iEta1 + 5;
0263 else if (point.eta() < -2.868 && point.eta() >= -2.871)
0264 return iEta1 + 4;
0265 else if (point.eta() < -2.871 && point.eta() >= -2.958)
0266 return iEta1 + 5;
0267 else if (point.eta() < -2.349)
0268 return iEta1 + 6;
0269 else
0270 return iEta1;
0271 }
0272
0273
0274 int HDetIdAssociator::iPhi(const GlobalPoint& point) {
0275 double pi = 4 * atan(1.);
0276 int iPhi1 = int((double(point.phi()) + pi) / (2 * pi) * nPhi_);
0277 return iPhi1;
0278 }
0279
0280
0281 void HDetIdAssociator::buildMap() {
0282
0283 check_setup();
0284 LogTrace("HDetIdAssociator") << "building map"
0285 << "\n";
0286 if (theMap_)
0287 delete theMap_;
0288 theMap_ = new std::vector<std::vector<std::set<DetId> > >(nEta_, std::vector<std::set<DetId> >(nPhi_));
0289 int numberOfDetIdsOutsideEtaRange = 0;
0290 int numberOfDetIdsActive = 0;
0291 std::set<DetId> validIds = getASetOfValidDetIds();
0292 for (std::set<DetId>::const_iterator id_itr = validIds.begin(); id_itr != validIds.end(); id_itr++) {
0293
0294 GlobalPoint point = getPosition(*id_itr);
0295
0296 if (point.eta() == 0)
0297 continue;
0298
0299 int ieta = iEta(point);
0300 int iphi = iPhi(point);
0301 int etaMax(-1);
0302 int etaMin(nEta_);
0303 int phiMax(-1);
0304 int phiMin(nPhi_);
0305 if (iphi >= nPhi_)
0306 iphi = iphi % nPhi_;
0307 assert(iphi >= 0);
0308 if (etaMin > ieta)
0309 etaMin = ieta;
0310 if (etaMax < ieta)
0311 etaMax = ieta;
0312 if (phiMin > iphi)
0313 phiMin = iphi;
0314 if (phiMax < iphi)
0315 phiMax = iphi;
0316
0317 if ((ieta > 54 || ieta < 15) && iphi % 2 == 0)
0318 phiMax++;
0319 if ((ieta > 54 || ieta < 15) && iphi % 2 == 1)
0320 phiMin--;
0321
0322 if (etaMax < 0 || phiMax < 0 || etaMin >= nEta_ || phiMin >= nPhi_) {
0323 LogTrace("HDetIdAssociator") << "Out of range: DetId:" << id_itr->rawId() << "\n\teta (min,max): " << etaMin
0324 << "," << etaMax << "\n\tphi (min,max): " << phiMin << "," << phiMax
0325 << "\nTower id: " << id_itr->rawId() << "\n";
0326 numberOfDetIdsOutsideEtaRange++;
0327 continue;
0328 }
0329
0330 if (phiMax - phiMin > phiMin + nPhi_ - phiMax) {
0331 phiMin += nPhi_;
0332 std::swap(phiMin, phiMax);
0333 }
0334 for (int ieta = etaMin; ieta <= etaMax; ieta++)
0335 for (int iphi = phiMin; iphi <= phiMax; iphi++)
0336 (*theMap_)[ieta][iphi % nPhi_].insert(*id_itr);
0337 numberOfDetIdsActive++;
0338 }
0339 LogTrace("HDetIdAssociator") << "Number of elements outside the allowed range ( |eta|>" << nEta_ / 2 * etaBinSize_
0340 << "): " << numberOfDetIdsOutsideEtaRange << "\n";
0341 LogTrace("HDetIdAssociator") << "Number of active DetId's mapped: " << numberOfDetIdsActive << "\n";
0342 }
0343
0344
0345 std::set<DetId> HDetIdAssociator::getDetIdsInACone(const std::set<DetId>& inset,
0346 const std::vector<GlobalPoint>& trajectory,
0347 const double dR) {
0348
0349 check_setup();
0350 std::set<DetId> outset;
0351
0352 if (dR >= 0) {
0353 for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
0354 for (std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end();
0355 point_iter++)
0356 if (nearElement(*point_iter, *id_iter, dR))
0357 outset.insert(*id_iter);
0358 } else {
0359 if (inset.size() != 1)
0360 return outset;
0361 std::set<DetId>::const_iterator id_inp = inset.begin();
0362 int ieta;
0363 int iphi;
0364 GlobalPoint point = getPosition(*id_inp);
0365 ieta = iEta(point);
0366 iphi = iPhi(point);
0367 for (int i = ieta - 1; i <= ieta + 1; i++) {
0368 for (int j = iphi - 1; j <= iphi + 1; j++) {
0369
0370 if (i < 0 || i >= nEta_)
0371 continue;
0372 int j2fill = j % nPhi_;
0373 if (j2fill < 0)
0374 j2fill += nPhi_;
0375 if ((*theMap_)[i][j2fill].empty())
0376 continue;
0377 outset.insert((*theMap_)[i][j2fill].begin(), (*theMap_)[i][j2fill].end());
0378 }
0379 }
0380 }
0381
0382
0383
0384
0385
0386
0387
0388
0389
0390 return outset;
0391 }
0392
0393
0394 std::set<DetId> HDetIdAssociator::getCrossedDetIds(const std::set<DetId>& inset,
0395 const std::vector<GlobalPoint>& trajectory) {
0396 check_setup();
0397 std::set<DetId> outset;
0398 for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++)
0399 for (std::vector<GlobalPoint>::const_iterator point_iter = trajectory.begin(); point_iter != trajectory.end();
0400 point_iter++)
0401 if (insideElement(*point_iter, *id_iter))
0402 outset.insert(*id_iter);
0403 return outset;
0404 }
0405
0406
0407 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
0408 edm::Handle<CaloTowerCollection> caloTowers) {
0409
0410 check_setup();
0411 std::set<DetId> outset;
0412 std::set<DetId>::const_iterator id_max = inset.begin();
0413 double Ehadmax = 0;
0414
0415 for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
0416 DetId id(*id_iter);
0417
0418
0419
0420 CaloTowerCollection::const_iterator tower = (*caloTowers).find(id);
0421 if (tower != (*caloTowers).end() && tower->hadEnergy() > Ehadmax) {
0422 id_max = id_iter;
0423 Ehadmax = tower->hadEnergy();
0424 }
0425 }
0426
0427 if (Ehadmax > 0)
0428 outset.insert(*id_max);
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438 return outset;
0439 }
0440
0441
0442 std::set<DetId> HDetIdAssociator::getMaxEDetId(const std::set<DetId>& inset,
0443 edm::Handle<HBHERecHitCollection> recHits) {
0444
0445 check_setup();
0446 std::set<DetId> outset;
0447 std::set<DetId>::const_iterator id_max = inset.begin();
0448 double Ehadmax = 0;
0449
0450 for (std::set<DetId>::const_iterator id_iter = inset.begin(); id_iter != inset.end(); id_iter++) {
0451 DetId id(*id_iter);
0452
0453
0454
0455 HBHERecHitCollection::const_iterator hit = (*recHits).find(id);
0456 if (hit != (*recHits).end() && hit->energy() > Ehadmax) {
0457 id_max = id_iter;
0458 Ehadmax = hit->energy();
0459 }
0460 }
0461
0462 if (Ehadmax > 0)
0463 outset.insert(*id_max);
0464
0465
0466
0467
0468
0469
0470
0471
0472
0473 return outset;
0474 }