File indexing completed on 2024-09-07 04:37:54
0001 #include "RecoPPS/Local/interface/RPixPlaneCombinatoryTracking.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <memory>
0007
0008 #include "DataFormats/Math/interface/Error.h"
0009 #include "DataFormats/Math/interface/AlgebraicROOTObjects.h"
0010 #include "CondFormats/PPSObjects/interface/CTPPSPixelIndices.h"
0011 #include "TMath.h"
0012
0013 #include "DataFormats/CTPPSReco/interface/CTPPSPixelLocalTrackRecoInfo.h"
0014
0015
0016
0017 RPixPlaneCombinatoryTracking::RPixPlaneCombinatoryTracking(edm::ParameterSet const ¶meterSet)
0018 : RPixDetTrackFinder(parameterSet) {
0019 trackMinNumberOfPoints_ = parameterSet.getParameter<uint>("trackMinNumberOfPoints");
0020 verbosity_ = parameterSet.getUntrackedParameter<int>("verbosity");
0021 maximumChi2OverNDF_ = parameterSet.getParameter<double>("maximumChi2OverNDF");
0022 maximumXLocalDistanceFromTrack_ = parameterSet.getParameter<double>("maximumXLocalDistanceFromTrack");
0023 maximumYLocalDistanceFromTrack_ = parameterSet.getParameter<double>("maximumYLocalDistanceFromTrack");
0024 numberOfPlanesPerPot_ = parameterSet.getParameter<int>("numberOfPlanesPerPot");
0025
0026 if (trackMinNumberOfPoints_ < 3) {
0027 throw cms::Exception("RPixPlaneCombinatoryTracking")
0028 << "Minimum number of planes required for tracking is 3, "
0029 << "tracking is not possible with " << trackMinNumberOfPoints_ << " hits";
0030 }
0031 }
0032
0033
0034
0035 RPixPlaneCombinatoryTracking::~RPixPlaneCombinatoryTracking() { possiblePlaneCombinations_.clear(); }
0036
0037
0038
0039 void RPixPlaneCombinatoryTracking::initialize() {
0040 uint32_t numberOfCombinations =
0041 factorial(numberOfPlanesPerPot_) /
0042 (factorial(numberOfPlanesPerPot_ - trackMinNumberOfPoints_) * factorial(trackMinNumberOfPoints_));
0043 if (verbosity_ >= 2)
0044 edm::LogInfo("RPixPlaneCombinatoryTracking") << "Number of combinations = " << numberOfCombinations;
0045 possiblePlaneCombinations_.reserve(numberOfCombinations);
0046
0047 getPlaneCombinations(listOfAllPlanes_, trackMinNumberOfPoints_, possiblePlaneCombinations_);
0048
0049 if (verbosity_ >= 2) {
0050 for (const auto &vec : possiblePlaneCombinations_) {
0051 for (const auto &num : vec) {
0052 edm::LogInfo("RPixPlaneCombinatoryTracking") << num << " - ";
0053 }
0054 edm::LogInfo("RPixPlaneCombinatoryTracking");
0055 }
0056 }
0057 }
0058
0059
0060
0061
0062 void RPixPlaneCombinatoryTracking::getPlaneCombinations(const std::vector<uint32_t> &inputPlaneList,
0063 uint32_t numberToExtract,
0064 PlaneCombinations &planeCombinations) const {
0065 uint32_t numberOfPlanes = inputPlaneList.size();
0066 std::string bitmask(numberToExtract, 1);
0067 bitmask.resize(numberOfPlanes, 0);
0068 planeCombinations.clear();
0069
0070
0071 do {
0072 planeCombinations.emplace_back();
0073 for (uint32_t i = 0; i < numberOfPlanes; ++i) {
0074 if (bitmask[i])
0075 planeCombinations.back().push_back(inputPlaneList.at(i));
0076 }
0077 } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
0078
0079 return;
0080 }
0081
0082
0083
0084
0085
0086
0087
0088 void RPixPlaneCombinatoryTracking::getHitCombinations(const std::map<CTPPSPixelDetId, PointInPlaneList> &mapOfAllHits,
0089 std::map<CTPPSPixelDetId, PointInPlaneList>::iterator mapIterator,
0090 HitReferences tmpHitPlaneMap,
0091 const PointInPlaneList &tmpHitVector,
0092 PointAndReferenceMap &outputMap) {
0093
0094 if (mapIterator == mapOfAllHits.end()) {
0095 outputMap[tmpHitPlaneMap] = tmpHitVector;
0096 return;
0097 }
0098 for (size_t i = 0; i < mapIterator->second.size(); i++) {
0099 HitReferences newHitPlaneMap = tmpHitPlaneMap;
0100 newHitPlaneMap[mapIterator->first] = i;
0101 PointInPlaneList newVector = tmpHitVector;
0102 newVector.push_back(mapIterator->second[i]);
0103 std::map<CTPPSPixelDetId, PointInPlaneList>::iterator tmpMapIterator = mapIterator;
0104 getHitCombinations(mapOfAllHits, ++tmpMapIterator, newHitPlaneMap, newVector, outputMap);
0105 }
0106 }
0107
0108
0109
0110 RPixPlaneCombinatoryTracking::PointAndReferenceMap RPixPlaneCombinatoryTracking::produceAllHitCombination(
0111 PlaneCombinations inputPlaneCombination) {
0112 PointAndReferenceMap mapOfAllPoints;
0113 CTPPSPixelDetId tmpRpId = romanPotId_;
0114
0115 if (verbosity_ >= 2)
0116 edm::LogInfo("RPixPlaneCombinatoryTracking") << "Searching for all combinations...";
0117
0118 for (const auto &planeCombination : inputPlaneCombination) {
0119 std::map<CTPPSPixelDetId, PointInPlaneList> selectedCombinationHitOnPlane;
0120 bool allPlaneAsHits = true;
0121
0122
0123
0124
0125 for (const auto &plane : planeCombination) {
0126 tmpRpId.setPlane(plane);
0127 CTPPSPixelDetId planeDetId = tmpRpId;
0128 if (hitMap_->find(planeDetId) == hitMap_->end()) {
0129 if (verbosity_ >= 2)
0130 edm::LogInfo("RPixPlaneCombinatoryTracking")
0131 << "No data on arm " << planeDetId.arm() << " station " << planeDetId.station() << " rp "
0132 << planeDetId.rp() << " plane " << planeDetId.plane();
0133 allPlaneAsHits = false;
0134 break;
0135 }
0136 if (selectedCombinationHitOnPlane.find(planeDetId) != selectedCombinationHitOnPlane.end()) {
0137 throw cms::Exception("RPixPlaneCombinatoryTracking")
0138 << "selectedCombinationHitOnPlane contains already detId " << planeDetId
0139 << "Error in the algorithm which created all the possible plane combinations";
0140 }
0141 selectedCombinationHitOnPlane[planeDetId] = (*hitMap_)[planeDetId];
0142 }
0143 if (!allPlaneAsHits)
0144 continue;
0145
0146
0147 auto mapIterator = selectedCombinationHitOnPlane.begin();
0148 HitReferences tmpHitPlaneMap;
0149 PointInPlaneList tmpHitVector;
0150 getHitCombinations(selectedCombinationHitOnPlane, mapIterator, tmpHitPlaneMap, tmpHitVector, mapOfAllPoints);
0151 if (verbosity_ >= 2)
0152 edm::LogInfo("RPixPlaneCombinatoryTracking") << "Number of possible tracks " << mapOfAllPoints.size();
0153
0154 }
0155
0156 return mapOfAllPoints;
0157 }
0158
0159
0160
0161 void RPixPlaneCombinatoryTracking::findTracks(int run) {
0162
0163
0164
0165
0166
0167
0168 unsigned int hitNum = 0;
0169 for (const auto &plane : *hitMap_) {
0170 hitNum += plane.second.size();
0171 }
0172
0173 if (hitMap_->size() == 2 && hitNum == 2) {
0174
0175 GlobalPoint hit[2];
0176 PointInPlaneList pIPL;
0177
0178 unsigned int i = 0;
0179 for (const auto &plane : *hitMap_) {
0180 if (plane.second.size() > 1)
0181 break;
0182 GlobalPoint gP(plane.second[0].globalPoint.x(), plane.second[0].globalPoint.y(), plane.second[0].globalPoint.z());
0183 hit[i] = gP;
0184 i++;
0185 pIPL.push_back(plane.second[0]);
0186 }
0187
0188
0189 double tx = (hit[0].x() - hit[1].x()) / (hit[0].z() - hit[1].z());
0190 double ty = (hit[0].y() - hit[1].y()) / (hit[0].z() - hit[1].z());
0191 double xat0 = (hit[1].x() * hit[0].z() - hit[0].x() * hit[1].z()) / (hit[0].z() - hit[1].z());
0192 double yat0 = (hit[1].y() * hit[0].z() - hit[0].y() * hit[1].z()) / (hit[0].z() - hit[1].z());
0193 double z0 = geometry_->rpTranslation(romanPotId_).z();
0194 double xatz0 = xat0 + tx * z0;
0195 double yatz0 = yat0 + ty * z0;
0196
0197 math::Vector<4>::type parameterVector{xatz0, yatz0, tx, ty};
0198 ROOT::Math::SVector<double, 10> v(0.01, 0.0, 0.01, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0, 0.01);
0199 math::Error<4>::type covarianceMatrix(v);
0200
0201 CTPPSPixelLocalTrack track(z0, parameterVector, covarianceMatrix, 0);
0202
0203
0204 for (const auto &hhit : pIPL) {
0205 GlobalPoint pOD(hhit.globalPoint.x(), hhit.globalPoint.y(), hhit.globalPoint.z());
0206 LocalPoint res(0, 0);
0207 LocalPoint pulls(0, 0);
0208 CTPPSPixelFittedRecHit usedRecHit(hhit.recHit, pOD, res, pulls);
0209 usedRecHit.setIsRealHit(true);
0210 track.addHit(hhit.detId, usedRecHit);
0211 }
0212
0213
0214 localTrackVector_.push_back(track);
0215 }
0216
0217
0218 while (hitMap_->size() >= trackMinNumberOfPoints_) {
0219 if (verbosity_ >= 1)
0220 edm::LogInfo("RPixPlaneCombinatoryTracking") << "Number of plane with hits " << hitMap_->size();
0221 if (verbosity_ >= 2)
0222 for (const auto &plane : *hitMap_)
0223 edm::LogInfo("RPixPlaneCombinatoryTracking")
0224 << "\tarm " << plane.first.arm() << " station " << plane.first.station() << " rp " << plane.first.rp()
0225 << " plane " << plane.first.plane() << " : " << plane.second.size();
0226
0227
0228
0229 PointAndReferenceMap mapOfAllMinRequiredPoint;
0230
0231 mapOfAllMinRequiredPoint = produceAllHitCombination(possiblePlaneCombinations_);
0232
0233
0234 double theMinChiSquaredOverNDF = maximumChi2OverNDF_ + 1.;
0235 HitReferences pointMapWithMinChiSquared;
0236 PointInPlaneList pointsWithMinChiSquared;
0237 CTPPSPixelLocalTrack bestTrack;
0238
0239 if (verbosity_ >= 2)
0240 edm::LogInfo("RPixPlaneCombinatoryTracking")
0241 << "Number of combinations of trackMinNumberOfPoints_ planes " << mapOfAllMinRequiredPoint.size();
0242 for (const auto &pointsAndRef : mapOfAllMinRequiredPoint) {
0243 CTPPSPixelLocalTrack tmpTrack = fitTrack(pointsAndRef.second);
0244 double tmpChiSquaredOverNDF = tmpTrack.chiSquaredOverNDF();
0245 if (verbosity_ >= 2)
0246 edm::LogInfo("RPixPlaneCombinatoryTracking") << "ChiSquare of the present track " << tmpChiSquaredOverNDF;
0247 if (!tmpTrack.isValid() || tmpChiSquaredOverNDF > maximumChi2OverNDF_ || tmpChiSquaredOverNDF == 0.)
0248 continue;
0249 if (tmpChiSquaredOverNDF < theMinChiSquaredOverNDF) {
0250 theMinChiSquaredOverNDF = tmpChiSquaredOverNDF;
0251 pointMapWithMinChiSquared = pointsAndRef.first;
0252 pointsWithMinChiSquared = pointsAndRef.second;
0253 bestTrack = tmpTrack;
0254 }
0255 }
0256
0257
0258
0259 if (theMinChiSquaredOverNDF > maximumChi2OverNDF_)
0260 break;
0261
0262
0263 std::vector<uint32_t> listOfExcludedPlanes;
0264 for (const auto &plane : listOfAllPlanes_) {
0265 CTPPSPixelDetId tmpRpId = romanPotId_;
0266 tmpRpId.setPlane(plane);
0267 CTPPSPixelDetId planeDetId = tmpRpId;
0268 if (pointMapWithMinChiSquared.find(planeDetId) == pointMapWithMinChiSquared.end())
0269 listOfExcludedPlanes.push_back(plane);
0270 }
0271
0272
0273
0274 PlaneCombinations planePointedHitListCombination;
0275 for (uint32_t i = 1; i <= listOfExcludedPlanes.size(); ++i) {
0276 PlaneCombinations tmpPlaneCombination;
0277 getPlaneCombinations(listOfExcludedPlanes, i, tmpPlaneCombination);
0278 for (const auto &combination : tmpPlaneCombination)
0279 planePointedHitListCombination.push_back(combination);
0280 }
0281
0282
0283 PointAndReferenceMap mapOfAllPointWithAtLeastBestFitSelected;
0284 PointAndReferenceMap mapOfPointCombinationToBeAdded;
0285 mapOfPointCombinationToBeAdded = produceAllHitCombination(planePointedHitListCombination);
0286
0287 for (const auto &element : mapOfPointCombinationToBeAdded) {
0288 HitReferences newPointMap = pointMapWithMinChiSquared;
0289 PointInPlaneList newPoints = pointsWithMinChiSquared;
0290 for (const auto &pointRef : element.first)
0291 newPointMap[pointRef.first] = pointRef.second;
0292 for (const auto &point : element.second)
0293 newPoints.push_back(point);
0294 mapOfAllPointWithAtLeastBestFitSelected[newPointMap] = newPoints;
0295 }
0296
0297
0298 if (verbosity_ >= 1)
0299 edm::LogInfo("RPixPlaneCombinatoryTracking")
0300 << "Minimum chiSquare over NDF for all the tracks " << theMinChiSquaredOverNDF;
0301
0302
0303
0304 std::vector<PointAndReferencePair> orderedVectorOfAllPointWithAtLeastBestFitSelected =
0305 orderCombinationsPerNumberOrPoints(mapOfAllPointWithAtLeastBestFitSelected);
0306 int currentNumberOfPlanes = 0;
0307 theMinChiSquaredOverNDF = maximumChi2OverNDF_ + 1.;
0308 bool foundTrackWithCurrentNumberOfPlanes = false;
0309 for (const auto &pointsAndRef : orderedVectorOfAllPointWithAtLeastBestFitSelected) {
0310 int tmpNumberOfPlanes = pointsAndRef.second.size();
0311
0312 if (foundTrackWithCurrentNumberOfPlanes && tmpNumberOfPlanes < currentNumberOfPlanes)
0313 break;
0314 CTPPSPixelLocalTrack tmpTrack = fitTrack(pointsAndRef.second);
0315 double tmpChiSquaredOverNDF = tmpTrack.chiSquaredOverNDF();
0316 if (!tmpTrack.isValid() || tmpChiSquaredOverNDF > maximumChi2OverNDF_ || tmpChiSquaredOverNDF == 0.)
0317 continue;
0318 if (tmpChiSquaredOverNDF < theMinChiSquaredOverNDF) {
0319 theMinChiSquaredOverNDF = tmpChiSquaredOverNDF;
0320 pointMapWithMinChiSquared = pointsAndRef.first;
0321 bestTrack = tmpTrack;
0322 currentNumberOfPlanes = tmpNumberOfPlanes;
0323 foundTrackWithCurrentNumberOfPlanes = true;
0324 }
0325 }
0326
0327 if (verbosity_ >= 1)
0328 edm::LogInfo("RPixPlaneCombinatoryTracking") << "The best track has " << bestTrack.ndf() / 2 + 2;
0329
0330 std::vector<uint32_t> listOfPlaneNotUsedForFit = listOfAllPlanes_;
0331
0332 for (const auto &hitToErase : pointMapWithMinChiSquared) {
0333 std::map<CTPPSPixelDetId, PointInPlaneList>::iterator hitMapElement = hitMap_->find(hitToErase.first);
0334 if (hitMapElement == hitMap_->end()) {
0335 throw cms::Exception("RPixPlaneCombinatoryTracking")
0336 << "The found tracks has hit belonging to a plane which does not have hits";
0337 }
0338 std::vector<uint32_t>::iterator planeIt;
0339 planeIt = std::find(listOfPlaneNotUsedForFit.begin(), listOfPlaneNotUsedForFit.end(), hitToErase.first.plane());
0340 listOfPlaneNotUsedForFit.erase(planeIt);
0341 hitMapElement->second.erase(hitMapElement->second.begin() + hitToErase.second);
0342
0343 if (hitMapElement->second.empty())
0344 hitMap_->erase(hitMapElement);
0345 }
0346
0347
0348
0349
0350
0351
0352
0353 for (const auto &plane : listOfPlaneNotUsedForFit) {
0354 CTPPSPixelDetId tmpPlaneId = romanPotId_;
0355 tmpPlaneId.setPlane(plane);
0356 std::unique_ptr<CTPPSPixelFittedRecHit> fittedRecHit(new CTPPSPixelFittedRecHit());
0357 GlobalPoint pointOnDet;
0358 calculatePointOnDetector(&bestTrack, tmpPlaneId, pointOnDet);
0359
0360 if (hitMap_->find(tmpPlaneId) != hitMap_->end()) {
0361
0362
0363 math::Vector<3>::type maxGlobalPointDistance(
0364 maximumXLocalDistanceFromTrack_, maximumYLocalDistanceFromTrack_, 0.);
0365
0366 DetGeomDesc::RotationMatrix theRotationMatrix = geometry_->sensor(tmpPlaneId)->rotation();
0367 AlgebraicMatrix33 tmpPlaneRotationMatrixMap;
0368 theRotationMatrix.GetComponents(tmpPlaneRotationMatrixMap(0, 0),
0369 tmpPlaneRotationMatrixMap(0, 1),
0370 tmpPlaneRotationMatrixMap(0, 2),
0371 tmpPlaneRotationMatrixMap(1, 0),
0372 tmpPlaneRotationMatrixMap(1, 1),
0373 tmpPlaneRotationMatrixMap(1, 2),
0374 tmpPlaneRotationMatrixMap(2, 0),
0375 tmpPlaneRotationMatrixMap(2, 1),
0376 tmpPlaneRotationMatrixMap(2, 2));
0377
0378 maxGlobalPointDistance = tmpPlaneRotationMatrixMap * maxGlobalPointDistance;
0379
0380 double maximumXdistance = maxGlobalPointDistance[0] * maxGlobalPointDistance[0];
0381 double maximumYdistance = maxGlobalPointDistance[1] * maxGlobalPointDistance[1];
0382
0383 double minimumDistance = 1. + maximumXdistance + maximumYdistance;
0384 for (const auto &hit : (*hitMap_)[tmpPlaneId]) {
0385 double xResidual = hit.globalPoint.x() - pointOnDet.x();
0386 double yResidual = hit.globalPoint.y() - pointOnDet.y();
0387 double xDistance = xResidual * xResidual;
0388 double yDistance = yResidual * yResidual;
0389 double distance = xDistance + yDistance;
0390 if (xDistance < maximumXdistance && yDistance < maximumYdistance && distance < minimumDistance) {
0391 LocalPoint residuals(xResidual, yResidual, 0.);
0392 math::Error<3>::type globalError = hit.globalError;
0393 LocalPoint pulls(xResidual / std::sqrt(globalError[0][0]), yResidual / std::sqrt(globalError[1][1]), 0.);
0394 fittedRecHit = std::make_unique<CTPPSPixelFittedRecHit>(hit.recHit, pointOnDet, residuals, pulls);
0395 fittedRecHit->setIsRealHit(true);
0396 }
0397 }
0398 } else {
0399 LocalPoint fakePoint;
0400 LocalError fakeError;
0401 CTPPSPixelRecHit fakeRecHit(fakePoint, fakeError);
0402 fittedRecHit = std::make_unique<CTPPSPixelFittedRecHit>(fakeRecHit, pointOnDet, fakePoint, fakePoint);
0403 }
0404
0405 bestTrack.addHit(tmpPlaneId, *fittedRecHit);
0406 }
0407
0408 localTrackVector_.push_back(bestTrack);
0409
0410 int pointForTracking = 0;
0411 int pointOnTrack = 0;
0412
0413 if (verbosity_ >= 1) {
0414 for (const auto &planeHits : bestTrack.hits()) {
0415 for (const auto &fittedhit : planeHits) {
0416 if (fittedhit.isUsedForFit())
0417 ++pointForTracking;
0418 if (fittedhit.isRealHit())
0419 ++pointOnTrack;
0420 }
0421 }
0422 edm::LogInfo("RPixPlaneCombinatoryTracking")
0423 << "Best track has " << pointForTracking << " points used for the fit and " << pointOnTrack
0424 << " points belonging to the track\n";
0425 }
0426
0427 }
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440 const CTPPSPixelDetId rpId_arm0_st2 = CTPPSPixelDetId(0, 2, 3);
0441 const CTPPSPixelDetId rpId_arm1_st2 = CTPPSPixelDetId(1, 2, 3);
0442 static const std::map<unsigned int, std::map<CTPPSPixelDetId, std::vector<bool> > > isPlaneShifted = {
0443 {0,
0444 {
0445 {rpId_arm0_st2, {false, false, false, false, false, false}},
0446 {rpId_arm1_st2, {false, false, false, false, false, false}}
0447 }},
0448 {300802,
0449 {
0450 {rpId_arm0_st2, {true, false, true, true, false, false}},
0451 {rpId_arm1_st2, {false, false, false, false, false, false}}
0452 }},
0453 {303338,
0454 {
0455 {rpId_arm0_st2, {false, false, false, false, false, false}},
0456 {rpId_arm1_st2, {false, false, false, false, false, false}}
0457 }},
0458 {305169,
0459 {
0460 {rpId_arm0_st2, {false, true, false, true, false, true}},
0461 {rpId_arm1_st2, {false, false, false, false, false, false}}
0462 }},
0463 {305965,
0464 {
0465 {rpId_arm0_st2, {false, true, false, true, false, true}},
0466 {rpId_arm1_st2, {false, false, true, false, true, true}}
0467 }},
0468 {307083,
0469 {
0470 {rpId_arm0_st2, {false, false, false, false, false, false}},
0471 {rpId_arm1_st2, {false, false, false, false, false, false}}
0472 }}};
0473 const auto &shiftStatusInitialRun = std::prev(isPlaneShifted.upper_bound(run));
0474 unsigned short shiftedROC = 10;
0475 CTPPSPixelIndices pixelIndices;
0476
0477
0478 if (romanPotId_.arm() == 0)
0479 shiftedROC = 0;
0480 if (romanPotId_.arm() == 1)
0481 shiftedROC = 5;
0482
0483
0484 for (auto &track : localTrackVector_) {
0485 if (romanPotId_ != rpId_arm0_st2 && romanPotId_ != rpId_arm1_st2) {
0486 track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::notShiftedRun);
0487 if (verbosity_ >= 2)
0488 edm::LogInfo("RPixPlaneCombinatoryTracking") << "Analyzing run: " << run << "\nTrack belongs to Arm "
0489 << romanPotId_.arm() << " Station " << romanPotId_.station();
0490
0491 continue;
0492 }
0493 unsigned short bxShiftedPlanesUsed = 0;
0494 unsigned short bxNonShiftedPlanesUsed = 0;
0495 unsigned short hitInShiftedROC = 0;
0496
0497 auto const &fittedHits = track.hits();
0498 auto const &planeFlags = (shiftStatusInitialRun->second).at(romanPotId_);
0499
0500 for (const auto &planeHits : fittedHits) {
0501 unsigned short plane = CTPPSPixelDetId(planeHits.detId()).plane();
0502 for (const auto &hit : planeHits) {
0503 if (hit.isUsedForFit()) {
0504 if (pixelIndices.getROCId(hit.minPixelCol(), hit.minPixelRow()) == shiftedROC)
0505 hitInShiftedROC++;
0506 if (planeFlags.at(plane))
0507 bxShiftedPlanesUsed++;
0508 else if (planeFlags != std::vector<bool>(6, false))
0509 bxNonShiftedPlanesUsed++;
0510 }
0511 }
0512 }
0513
0514
0515 track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
0516 invalid);
0517 if (hitInShiftedROC < 3)
0518 track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::notShiftedRun);
0519 else {
0520 if (bxShiftedPlanesUsed == 0 && bxNonShiftedPlanesUsed == 0)
0521 track.setRecoInfo(
0522 CTPPSpixelLocalTrackReconstructionInfo::notShiftedRun);
0523 if (bxShiftedPlanesUsed == 3 && bxNonShiftedPlanesUsed == 0)
0524 track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
0525 allShiftedPlanes);
0526 if (bxShiftedPlanesUsed == 0 && bxNonShiftedPlanesUsed == 3)
0527
0528 track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::noShiftedPlanes);
0529 if (bxShiftedPlanesUsed > 0 && bxNonShiftedPlanesUsed > 0)
0530 track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
0531 mixedPlanes);
0532 }
0533 if (bxShiftedPlanesUsed + bxNonShiftedPlanesUsed > 6) {
0534 throw cms::Exception("RPixPlaneCombinatoryTracking") << "Error in RPixPlaneCombinatoryTracking::findTracks -> "
0535 << "More than six points found for a track, skipping.";
0536 continue;
0537 }
0538 if (track.recoInfo() == CTPPSpixelLocalTrackReconstructionInfo::invalid) {
0539 throw cms::Exception("RPixPlaneCombinatoryTracking") << "Error in RPixPlaneCombinatoryTracking::findTracks -> "
0540 << "recoInfo has not been set properly.";
0541 }
0542
0543 if (verbosity_ >= 2) {
0544 edm::LogInfo("RPixPlaneCombinatoryTracking")
0545 << " Track belongs to Arm " << romanPotId_.arm() << " Station " << romanPotId_.station()
0546 << "\nFirst run with this bx-shift configuration: " << shiftStatusInitialRun->first
0547 << "\nTrack reconstructed with: " << bxShiftedPlanesUsed << " bx-shifted planes, " << bxNonShiftedPlanesUsed
0548 << " non-bx-shifted planes, " << hitInShiftedROC << " hits in the bx-shifted ROC"
0549 << "\nrecoInfo = " << (unsigned short)track.recoInfo();
0550 if (planeFlags != std::vector<bool>(6, false))
0551 edm::LogInfo("RPixPlaneCombinatoryTracking") << "The shifted ROC is ROC" << shiftedROC;
0552 }
0553 }
0554 return;
0555 }
0556
0557
0558
0559 CTPPSPixelLocalTrack RPixPlaneCombinatoryTracking::fitTrack(PointInPlaneList pointList) {
0560 uint32_t const numberOfPlanes = 6;
0561 math::Vector<2 * numberOfPlanes>::type xyCoordinates;
0562 math::Error<2 * numberOfPlanes>::type varianceMatrix;
0563 math::Matrix<2 * numberOfPlanes, 4>::type zMatrix;
0564
0565
0566 for (uint32_t iHit = 0; iHit < numberOfPlanes; iHit++) {
0567 if (iHit < pointList.size()) {
0568 const auto &globalPoint = pointList[iHit].globalPoint;
0569 xyCoordinates[2 * iHit] = globalPoint.x();
0570 xyCoordinates[2 * iHit + 1] = globalPoint.y();
0571 zMatrix(2 * iHit, 0) = 1.;
0572 zMatrix(2 * iHit, 2) = globalPoint.z() - z0_;
0573 zMatrix(2 * iHit + 1, 1) = 1.;
0574 zMatrix(2 * iHit + 1, 3) = globalPoint.z() - z0_;
0575
0576 AlgebraicMatrix33 globalError = pointList[iHit].globalError;
0577 varianceMatrix(2 * iHit, 2 * iHit) = globalError(0, 0);
0578 varianceMatrix(2 * iHit, 2 * iHit + 1) = globalError(0, 1);
0579 varianceMatrix(2 * iHit + 1, 2 * iHit) = globalError(1, 0);
0580 varianceMatrix(2 * iHit + 1, 2 * iHit + 1) = globalError(1, 1);
0581 } else {
0582 varianceMatrix(2 * iHit, 2 * iHit) = 1.;
0583 varianceMatrix(2 * iHit + 1, 2 * iHit + 1) = 1.;
0584 }
0585 }
0586
0587
0588 if (!varianceMatrix.Invert()) {
0589 edm::LogError("RPixPlaneCombinatoryTracking") << "Error in RPixPlaneCombinatoryTracking::fitTrack -> "
0590 << "Point variance matrix is singular, skipping.";
0591 CTPPSPixelLocalTrack badTrack;
0592 badTrack.setValid(false);
0593 return badTrack;
0594 }
0595
0596 math::Error<4>::type covarianceMatrix = ROOT::Math::SimilarityT(zMatrix, varianceMatrix);
0597
0598
0599 if (!covarianceMatrix.Invert()) {
0600 edm::LogError("RPixPlaneCombinatoryTracking") << "Error in RPixPlaneCombinatoryTracking::fitTrack -> "
0601 << "Parameter covariance matrix is singular, skipping.";
0602 CTPPSPixelLocalTrack badTrack;
0603 badTrack.setValid(false);
0604 return badTrack;
0605 }
0606
0607
0608 math::Vector<4>::type zMatrixTransposeTimesVarianceMatrixTimesXyCoordinates =
0609 ROOT::Math::Transpose(zMatrix) * varianceMatrix * xyCoordinates;
0610 math::Vector<4>::type parameterVector = covarianceMatrix * zMatrixTransposeTimesVarianceMatrixTimesXyCoordinates;
0611 math::Vector<2 * numberOfPlanes>::type xyCoordinatesMinusZmatrixTimesParameters =
0612 xyCoordinates - (zMatrix * parameterVector);
0613
0614 double chiSquare = ROOT::Math::Dot(xyCoordinatesMinusZmatrixTimesParameters,
0615 (varianceMatrix * xyCoordinatesMinusZmatrixTimesParameters));
0616
0617 CTPPSPixelLocalTrack goodTrack(z0_, parameterVector, covarianceMatrix, chiSquare);
0618 goodTrack.setValid(true);
0619
0620 for (const auto &hit : pointList) {
0621 const auto &globalPoint = hit.globalPoint;
0622 GlobalPoint pointOnDet;
0623 bool foundPoint = calculatePointOnDetector(&goodTrack, hit.detId, pointOnDet);
0624 if (!foundPoint) {
0625 CTPPSPixelLocalTrack badTrack;
0626 badTrack.setValid(false);
0627 return badTrack;
0628 }
0629 double xResidual = globalPoint.x() - pointOnDet.x();
0630 double yResidual = globalPoint.y() - pointOnDet.y();
0631 LocalPoint residuals(xResidual, yResidual);
0632
0633 math::Error<3>::type globalError(hit.globalError);
0634 LocalPoint pulls(xResidual / std::sqrt(globalError(0, 0)), yResidual / std::sqrt(globalError(0, 0)));
0635
0636 CTPPSPixelFittedRecHit fittedRecHit(hit.recHit, pointOnDet, residuals, pulls);
0637 fittedRecHit.setIsUsedForFit(true);
0638 goodTrack.addHit(hit.detId, fittedRecHit);
0639 }
0640
0641 return goodTrack;
0642 }
0643
0644
0645
0646
0647 bool RPixPlaneCombinatoryTracking::calculatePointOnDetector(CTPPSPixelLocalTrack *track,
0648 CTPPSPixelDetId planeId,
0649 GlobalPoint &planeLineIntercept) {
0650 double z0 = track->z0();
0651 CTPPSPixelLocalTrack::ParameterVector parameters = track->parameterVector();
0652
0653 math::Vector<3>::type pointOnLine(parameters[0], parameters[1], z0);
0654 GlobalVector tmpLineUnitVector = track->directionVector();
0655 math::Vector<3>::type lineUnitVector(tmpLineUnitVector.x(), tmpLineUnitVector.y(), tmpLineUnitVector.z());
0656
0657 const CTPPSGeometry::Vector tmpPointLocal(0., 0., 0.);
0658 const auto &tmpPointOnPlane = geometry_->localToGlobal(planeId, tmpPointLocal);
0659
0660 math::Vector<3>::type pointOnPlane(tmpPointOnPlane.x(), tmpPointOnPlane.y(), tmpPointOnPlane.z());
0661 math::Vector<3>::type planeUnitVector(0., 0., 1.);
0662
0663 DetGeomDesc::RotationMatrix theRotationMatrix = geometry_->sensor(planeId)->rotation();
0664 AlgebraicMatrix33 tmpPlaneRotationMatrixMap;
0665 theRotationMatrix.GetComponents(tmpPlaneRotationMatrixMap(0, 0),
0666 tmpPlaneRotationMatrixMap(0, 1),
0667 tmpPlaneRotationMatrixMap(0, 2),
0668 tmpPlaneRotationMatrixMap(1, 0),
0669 tmpPlaneRotationMatrixMap(1, 1),
0670 tmpPlaneRotationMatrixMap(1, 2),
0671 tmpPlaneRotationMatrixMap(2, 0),
0672 tmpPlaneRotationMatrixMap(2, 1),
0673 tmpPlaneRotationMatrixMap(2, 2));
0674
0675 planeUnitVector = tmpPlaneRotationMatrixMap * planeUnitVector;
0676
0677 double denominator = ROOT::Math::Dot(lineUnitVector, planeUnitVector);
0678 if (denominator == 0) {
0679 edm::LogError("RPixPlaneCombinatoryTracking")
0680 << "Error in RPixPlaneCombinatoryTracking::calculatePointOnDetector -> "
0681 << "Fitted line and plane are parallel. Removing this track";
0682 return false;
0683 }
0684
0685 double distanceFromLinePoint = ROOT::Math::Dot((pointOnPlane - pointOnLine), planeUnitVector) / denominator;
0686
0687 math::Vector<3>::type tmpPlaneLineIntercept = distanceFromLinePoint * lineUnitVector + pointOnLine;
0688 planeLineIntercept = GlobalPoint(tmpPlaneLineIntercept[0], tmpPlaneLineIntercept[1], tmpPlaneLineIntercept[2]);
0689
0690 return true;
0691 }
0692
0693
0694
0695 std::vector<RPixPlaneCombinatoryTracking::PointAndReferencePair>
0696 RPixPlaneCombinatoryTracking::orderCombinationsPerNumberOrPoints(PointAndReferenceMap inputMap) {
0697 std::vector<PointAndReferencePair> sortedVector(inputMap.begin(), inputMap.end());
0698 std::sort(sortedVector.begin(), sortedVector.end(), functionForPlaneOrdering);
0699
0700 return sortedVector;
0701 }
0702