Local/src/RPixPlaneCombinatoryTracking.cc

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 &parameterSet)
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 //This function produces all the possible plane combinations extracting numberToExtract planes over numberOfPlanes planes
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);  // numberToExtract leading 1's
0067   bitmask.resize(numberOfPlanes, 0);        // numberOfPlanes-numberToExtract trailing 0's
0068   planeCombinations.clear();
0069
0070   // store the combination and permute bitmask
0071   do {
0072     planeCombinations.emplace_back();
0073     for (uint32_t i = 0; i < numberOfPlanes; ++i) {  // [0..numberOfPlanes-1] integers
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 //This function calls it self in order to get all the possible combinations of hits having a certain selected number of planes
0085 //This function avoids to write for loops in cascade which will not allow to define arbitrarily the minimum number of planes to use
0086 //The output is stored in a map containing the vector of points and as a key the map of the point forming this vector.
0087 //This allows to erase the points already used for the track fit
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   //At this point I selected one hit per plane
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_;  //in order to avoid to modify the data member
0114
0115   if (verbosity_ >= 2)
0116     edm::LogInfo("RPixPlaneCombinatoryTracking") << "Searching for all combinations...";
0117   //Loop on all the plane combinations
0118   for (const auto &planeCombination : inputPlaneCombination) {
0119     std::map<CTPPSPixelDetId, PointInPlaneList> selectedCombinationHitOnPlane;
0120     bool allPlaneAsHits = true;
0121
0122     //Loop on all the possible combinations
0123     //In this loop the selectedCombinationHitOnPlane is filled
0124     //and cases in which one of the selected plane is empty are skipped
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     //I add the all the hit combinations to the full list of plane combinations
0147     auto mapIterator = selectedCombinationHitOnPlane.begin();
0148     HitReferences tmpHitPlaneMap;   //empty map of plane id and hit number needed the getHitCombinations algorithm
0149     PointInPlaneList tmpHitVector;  //empty vector of hits needed for the getHitCombinations algorithm
0150     getHitCombinations(selectedCombinationHitOnPlane, mapIterator, tmpHitPlaneMap, tmpHitVector, mapOfAllPoints);
0151     if (verbosity_ >= 2)
0152       edm::LogInfo("RPixPlaneCombinatoryTracking") << "Number of possible tracks " << mapOfAllPoints.size();
0153
0154   }  //end of loops on all the combinations
0155
0156   return mapOfAllPoints;
0157 }
0158
0159 //------------------------------------------------------------------------------------------------//
0160
0161 void RPixPlaneCombinatoryTracking::findTracks(int run) {
0162   //The loop search for all the possible tracks starting from the one with the smallest chiSquare/NDF
0163   //The loop stops when the number of planes with recorded hits is less than the minimum number of planes required
0164   //or if the track with minimum chiSquare found has a chiSquare higher than the maximum required
0165
0166   // bad Pot patch -- beginning
0167   // check number of hits in road
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) {  // look for roads with 2 hits in 2 different planes
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;  // only 1 hit per plane per road allowed
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     // calculate 2 point track parameters (no need of fits)
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     // save used points into track
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     // save track in collection
0214     localTrackVector_.push_back(track);
0215   }
0216   // bad Pot patch -- end
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     //I create the map of all the possible combinations of a group of trackMinNumberOfPoints_ points
0228     //and the key keeps the reference of which planes and which hit numbers form the combination
0229     PointAndReferenceMap mapOfAllMinRequiredPoint;
0230     //I produce the map for all cominations of all hits with all trackMinNumberOfPoints_ plane combinations
0231     mapOfAllMinRequiredPoint = produceAllHitCombination(possiblePlaneCombinations_);
0232
0233     //Fit all the possible combinations with minimum number of planes required and find the track with minimum chi2
0234     double theMinChiSquaredOverNDF = maximumChi2OverNDF_ + 1.;  //in order to break the loop in case no track is found;
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;  //validity check
0249       if (tmpChiSquaredOverNDF < theMinChiSquaredOverNDF) {
0250         theMinChiSquaredOverNDF = tmpChiSquaredOverNDF;
0251         pointMapWithMinChiSquared = pointsAndRef.first;
0252         pointsWithMinChiSquared = pointsAndRef.second;
0253         bestTrack = tmpTrack;
0254       }
0255     }
0256
0257     //The loop on the fit of all tracks is done, the track with minimum chiSquared is found
0258     // and it is verified that it complies with the maximumChi2OverNDF_ requirement
0259     if (theMinChiSquaredOverNDF > maximumChi2OverNDF_)
0260       break;
0261
0262     //The list of planes not included in the minimum chiSquare track is produced.
0263     std::vector<uint32_t> listOfExcludedPlanes;
0264     for (const auto &plane : listOfAllPlanes_) {
0265       CTPPSPixelDetId tmpRpId = romanPotId_;  //in order to avoid to modify the data member
0266       tmpRpId.setPlane(plane);
0267       CTPPSPixelDetId planeDetId = tmpRpId;
0268       if (pointMapWithMinChiSquared.find(planeDetId) == pointMapWithMinChiSquared.end())
0269         listOfExcludedPlanes.push_back(plane);
0270     }
0271
0272     //I produce all the possible combinations of planes to be added to the track,
0273     //excluding the case of no other plane added since it has been already fitted.
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     //I produce all the possible combinations of points to be added to the track
0283     PointAndReferenceMap mapOfAllPointWithAtLeastBestFitSelected;
0284     PointAndReferenceMap mapOfPointCombinationToBeAdded;
0285     mapOfPointCombinationToBeAdded = produceAllHitCombination(planePointedHitListCombination);
0286     //The found hit combination is added to the hits selected by the best fit;
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;  //add the new point reference
0292       for (const auto &point : element.second)
0293         newPoints.push_back(point);
0294       mapOfAllPointWithAtLeastBestFitSelected[newPointMap] = newPoints;
0295     }
0296
0297     //I fit all the possible combination of the minimum plane best fit hits plus hits from the other planes
0298     if (verbosity_ >= 1)
0299       edm::LogInfo("RPixPlaneCombinatoryTracking")
0300           << "Minimum chiSquare over NDF for all the tracks " << theMinChiSquaredOverNDF;
0301
0302     // I look for the tracks with maximum number of points with a chiSquare over NDF smaller than maximumChi2OverNDF_
0303     // If more than one track fulfill the chiSquare requirement with the same number of points I choose the one with smaller chiSquare
0304     std::vector<PointAndReferencePair> orderedVectorOfAllPointWithAtLeastBestFitSelected =
0305         orderCombinationsPerNumberOrPoints(mapOfAllPointWithAtLeastBestFitSelected);
0306     int currentNumberOfPlanes = 0;
0307     theMinChiSquaredOverNDF = maximumChi2OverNDF_ + 1.;  //in order to break the loop in case no track is found;
0308     bool foundTrackWithCurrentNumberOfPlanes = false;
0309     for (const auto &pointsAndRef : orderedVectorOfAllPointWithAtLeastBestFitSelected) {
0310       int tmpNumberOfPlanes = pointsAndRef.second.size();
0311       // If a valid track has been already found with an higher number of planes the loop stops.
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;  //validity check
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     //remove the hits belonging to the tracks from the full list of hits
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       //if the plane at which the hit was erased is empty it is removed from the hit map
0343       if (hitMapElement->second.empty())
0344         hitMap_->erase(hitMapElement);
0345     }
0346
0347     //search for hit on the other planes which may belong to the same track
0348     //even if they did not contributed to the track
0349     //in case of multiple hit, the closest one to the track will be considered
0350     //If a hit is found these will not be erased from the list of all hits
0351     //If no hit is found, the point on the plane intersecting the track will be saved by a CTPPSPixelFittedRecHit
0352     //with the isRealHit_ flag set to false
0353     for (const auto &plane : listOfPlaneNotUsedForFit) {
0354       CTPPSPixelDetId tmpPlaneId = romanPotId_;  //in order to avoid to modify the data member
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         //I convert the hit search window defined in local coordinated into global
0362         //This avoids to convert the global plane-line intersection in order not to call the the geometry
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         //I avoid the Sqrt since it will not be saved
0380         double maximumXdistance = maxGlobalPointDistance[0] * maxGlobalPointDistance[0];
0381         double maximumYdistance = maxGlobalPointDistance[1] * maxGlobalPointDistance[1];
0382         // to be sure that the first min distance is from a real point
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   }  //close of the while loop on all the hits
0428
0429   // recoInfo_ calculation
0430   // Hardcoded shift periods:
0431   // Before run 300802: No shift
0432   // Starting from run 300802: Sec45 St2 Rp3 Pl 0,2,3 ROC 0 shifted.
0433   // Starting from run 303338: No shift.
0434   // Starting from run 305169: Sec45 St2 Rp3 Pl 1,3,5 ROC 0 shifted.
0435   // Starting from run 305965: Sec45 St2 Rp3 Pl 1,3,5 ROC 0 shifted & Sec56 St2 Rp3 Pl 2,4,5 ROC 5 shifted.
0436   // Starting from run 307083: No shift
0437
0438   // These variables hold the information of the runs when the detector was taking data in 3+3 Mode and which planes were bx-shifted
0439   // These values will never be changed and the 3+3 Mode will never be used again in the future
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}},  // Shift Period 0 Sec45
0446            {rpId_arm1_st2, {false, false, false, false, false, false}}   // Shift Period 1 Sec56
0447        }},
0448       {300802,
0449        {
0450            {rpId_arm0_st2, {true, false, true, true, false, false}},    // Shift Period 1 Sec45
0451            {rpId_arm1_st2, {false, false, false, false, false, false}}  // Shift Period 1 Sec56
0452        }},
0453       {303338,
0454        {
0455            {rpId_arm0_st2, {false, false, false, false, false, false}},  // Shift Period 2 Sec45
0456            {rpId_arm1_st2, {false, false, false, false, false, false}}   // Shift Period 2 Sec56
0457        }},
0458       {305169,
0459        {
0460            {rpId_arm0_st2, {false, true, false, true, false, true}},    // Shift Period 3 Sec45
0461            {rpId_arm1_st2, {false, false, false, false, false, false}}  // Shift Period 3 Sec56
0462        }},
0463       {305965,
0464        {
0465            {rpId_arm0_st2, {false, true, false, true, false, true}},  // Shift Period 4 Sec45
0466            {rpId_arm1_st2, {false, false, true, false, true, true}}   // Shift Period 4 Sec56
0467        }},
0468       {307083,
0469        {
0470            {rpId_arm0_st2, {false, false, false, false, false, false}},  // Shift Period 0 Sec45
0471            {rpId_arm1_st2, {false, false, false, false, false, false}}   // Shift Period 1 Sec56
0472        }}};                                                              // map< shiftedPeriod, map<DetID,shiftScheme> >
0473   const auto &shiftStatusInitialRun = std::prev(isPlaneShifted.upper_bound(run));
0474   unsigned short shiftedROC = 10;
0475   CTPPSPixelIndices pixelIndices;
0476
0477   // Selecting the shifted ROC
0478   if (romanPotId_.arm() == 0)
0479     shiftedROC = 0;
0480   if (romanPotId_.arm() == 1)
0481     shiftedROC = 5;
0482
0483   // Loop over found tracks to set recoInfo_
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++;  // Count how many hits are in the shifted ROC
0506           if (planeFlags.at(plane))
0507             bxShiftedPlanesUsed++;  // Count how many bx-shifted planes are used
0508           else if (planeFlags != std::vector<bool>(6, false))
0509             bxNonShiftedPlanesUsed++;  // Count how many non-bx-shifted planes are used, only if there are shifted planes
0510         }
0511       }
0512     }
0513
0514     // Set recoInfo_ value
0515     track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
0516                           invalid);  // Initially setting it as invalid. It has to match one of the following options.
0517     if (hitInShiftedROC < 3)
0518       track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::notShiftedRun);
0519     else {
0520       if (bxShiftedPlanesUsed == 0 && bxNonShiftedPlanesUsed == 0)
0521         track.setRecoInfo(
0522             CTPPSpixelLocalTrackReconstructionInfo::notShiftedRun);  // Default value for runs without bx-shift
0523       if (bxShiftedPlanesUsed == 3 && bxNonShiftedPlanesUsed == 0)
0524         track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
0525                               allShiftedPlanes);  // Track reconstructed in a shifted ROC, only with bx-shifted planes
0526       if (bxShiftedPlanesUsed == 0 && bxNonShiftedPlanesUsed == 3)
0527         // Track reconstructed in a shifted ROC, only with non-bx-shifted planes
0528         track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::noShiftedPlanes);
0529       if (bxShiftedPlanesUsed > 0 && bxNonShiftedPlanesUsed > 0)
0530         track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
0531                               mixedPlanes);  // Track reconstructed in a shifted ROC, with mixed planes
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   //The matrices and vector xyCoordinates, varianceMatrix and varianceMatrix are built from the points
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   //Get real point variance matrix
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   //To have the real parameter covariance matrix, covarianceMatrix needs to be inverted
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   // track parameters: (x0, y0, tx, ty); x = x0 + tx*(z-z0)
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 //The method calculates the hit pointed by the track on the detector plane
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 // The method sorts the possible point combinations in order to process before fits on the highest possible number of points
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 //------------------------------------------------------------------------------------------------//