#include <RPixPlaneCombinatoryTracking.h>

Inheritance diagram for RPixPlaneCombinatoryTracking:

Public Member Functions
void	findTracks (int run) override

void	initialize () override

	RPixPlaneCombinatoryTracking (edm::ParameterSet const &parameterSet)

	~RPixPlaneCombinatoryTracking () override

Public Member Functions inherited from RPixDetTrackFinder
void	clear ()

const std::vector< CTPPSPixelLocalTrack > &	getLocalTracks () const

	RPixDetTrackFinder (edm::ParameterSet const &parameterSet)

void	setGeometry (const CTPPSGeometry *geometry)

void	setHits (std::map< CTPPSPixelDetId, std::vector< RPixDetPatternFinder::PointInPlane > > *hitMap)

void	setListOfPlanes (std::vector< uint32_t > listOfAllPlanes)

void	setRomanPotId (CTPPSPixelDetId rpId)

void	setZ0 (double z0)

virtual	~RPixDetTrackFinder ()

Private Types
typedef std::map< CTPPSPixelDetId, size_t >	HitReferences

typedef std::vector< std::vector< uint32_t > >	PlaneCombinations

typedef std::map< HitReferences, PointInPlaneList >	PointAndReferenceMap

typedef std::pair< HitReferences, PointInPlaneList >	PointAndReferencePair

typedef std::vector< RPixDetPatternFinder::PointInPlane >	PointInPlaneList

Private Member Functions
bool	calculatePointOnDetector (CTPPSPixelLocalTrack *track, CTPPSPixelDetId planeId, GlobalPoint &planeLineIntercept)

uint32_t	factorial (uint32_t x) const

CTPPSPixelLocalTrack	fitTrack (PointInPlaneList pointList)

void	getHitCombinations (const std::map< CTPPSPixelDetId, PointInPlaneList > &mapOfAllHits, std::map< CTPPSPixelDetId, PointInPlaneList >::iterator mapIterator, HitReferences tmpHitPlaneMap, const PointInPlaneList &tmpHitVector, PointAndReferenceMap &outputMap)

void	getPlaneCombinations (const std::vector< uint32_t > &inputPlaneList, uint32_t numberToExtract, PlaneCombinations &planeCombinations) const

std::vector< PointAndReferencePair >	orderCombinationsPerNumberOrPoints (PointAndReferenceMap inputMap)

PointAndReferenceMap	produceAllHitCombination (PlaneCombinations inputPlaneCombination)

Static Private Member Functions
static bool	functionForPlaneOrdering (PointAndReferencePair a, PointAndReferencePair b)

Private Attributes
double	maximumChi2OverNDF_

double	maximumXLocalDistanceFromTrack_

double	maximumYLocalDistanceFromTrack_

PlaneCombinations	possiblePlaneCombinations_

uint32_t	trackMinNumberOfPoints_

int	verbosity_

Additional Inherited Members
Protected Attributes inherited from RPixDetTrackFinder
const CTPPSGeometry *	geometry_

std::map< CTPPSPixelDetId, std::vector< RPixDetPatternFinder::PointInPlane > > *	hitMap_

std::vector< uint32_t >	listOfAllPlanes_

std::vector< CTPPSPixelLocalTrack >	localTrackVector_

uint32_t	numberOfPlanesPerPot_

CTPPSPixelDetId	romanPotId_

double	z0_

Detailed Description

Definition at line 24 of file RPixPlaneCombinatoryTracking.h.

Member Typedef Documentation

◆ HitReferences

typedef std::map<CTPPSPixelDetId, size_t> RPixPlaneCombinatoryTracking::HitReferences

private

Definition at line 34 of file RPixPlaneCombinatoryTracking.h.

◆ PlaneCombinations

typedef std::vector<std::vector<uint32_t> > RPixPlaneCombinatoryTracking::PlaneCombinations

private

Definition at line 32 of file RPixPlaneCombinatoryTracking.h.

◆ PointAndReferenceMap

typedef std::map<HitReferences, PointInPlaneList> RPixPlaneCombinatoryTracking::PointAndReferenceMap

private

Definition at line 35 of file RPixPlaneCombinatoryTracking.h.

◆ PointAndReferencePair

typedef std::pair<HitReferences, PointInPlaneList> RPixPlaneCombinatoryTracking::PointAndReferencePair

private

Definition at line 36 of file RPixPlaneCombinatoryTracking.h.

◆ PointInPlaneList

typedef std::vector<RPixDetPatternFinder::PointInPlane> RPixPlaneCombinatoryTracking::PointInPlaneList

private

Definition at line 33 of file RPixPlaneCombinatoryTracking.h.

Constructor & Destructor Documentation

◆ RPixPlaneCombinatoryTracking()

RPixPlaneCombinatoryTracking::RPixPlaneCombinatoryTracking ( edm::ParameterSet const & parameterSet )

Definition at line 16 of file RPixPlaneCombinatoryTracking.cc.

     : RPixDetTrackFinder(parameterSet) {
   trackMinNumberOfPoints_ = parameterSet.getParameter<uint>("trackMinNumberOfPoints");
   verbosity_ = parameterSet.getUntrackedParameter<int>("verbosity");
   maximumChi2OverNDF_ = parameterSet.getParameter<double>("maximumChi2OverNDF");
   maximumXLocalDistanceFromTrack_ = parameterSet.getParameter<double>("maximumXLocalDistanceFromTrack");
   maximumYLocalDistanceFromTrack_ = parameterSet.getParameter<double>("maximumYLocalDistanceFromTrack");
   numberOfPlanesPerPot_ = parameterSet.getParameter<int>("numberOfPlanesPerPot");
  
   if (trackMinNumberOfPoints_ < 3) {
     throw cms::Exception("RPixPlaneCombinatoryTracking")
         << "Minimum number of planes required for tracking is 3, "
         << "tracking is not possible with " << trackMinNumberOfPoints_ << " hits";
   }
 }

References Exception, edm::ParameterSet::getParameter(), edm::ParameterSet::getUntrackedParameter(), maximumChi2OverNDF_, maximumXLocalDistanceFromTrack_, maximumYLocalDistanceFromTrack_, RPixDetTrackFinder::numberOfPlanesPerPot_, edm::parameterSet(), trackMinNumberOfPoints_, parallelization::uint, and verbosity_.

◆ ~RPixPlaneCombinatoryTracking()

RPixPlaneCombinatoryTracking::~RPixPlaneCombinatoryTracking ( )

override

Definition at line 34 of file RPixPlaneCombinatoryTracking.cc.

34 { possiblePlaneCombinations_.clear(); }

References possiblePlaneCombinations_.

Member Function Documentation

◆ calculatePointOnDetector()

bool RPixPlaneCombinatoryTracking::calculatePointOnDetector	(	CTPPSPixelLocalTrack *	track,
		CTPPSPixelDetId	planeId,
		GlobalPoint &	planeLineIntercept
	)

private

Definition at line 594 of file RPixPlaneCombinatoryTracking.cc.

                                                                                              {
   double z0 = track->z0();
   CTPPSPixelLocalTrack::ParameterVector parameters = track->parameterVector();
  
   math::Vector<3>::type pointOnLine(parameters[0], parameters[1], z0);
   GlobalVector tmpLineUnitVector = track->directionVector();
   math::Vector<3>::type lineUnitVector(tmpLineUnitVector.x(), tmpLineUnitVector.y(), tmpLineUnitVector.z());
  
   const CTPPSGeometry::Vector tmpPointLocal(0., 0., 0.);
   const auto &tmpPointOnPlane = geometry_->localToGlobal(planeId, tmpPointLocal);
  
   math::Vector<3>::type pointOnPlane(tmpPointOnPlane.x(), tmpPointOnPlane.y(), tmpPointOnPlane.z());
   math::Vector<3>::type planeUnitVector(0., 0., 1.);
  
   DetGeomDesc::RotationMatrix theRotationMatrix = geometry_->sensor(planeId)->rotation();
   AlgebraicMatrix33 tmpPlaneRotationMatrixMap;
   theRotationMatrix.GetComponents(tmpPlaneRotationMatrixMap(0, 0),
                                   tmpPlaneRotationMatrixMap(0, 1),
                                   tmpPlaneRotationMatrixMap(0, 2),
                                   tmpPlaneRotationMatrixMap(1, 0),
                                   tmpPlaneRotationMatrixMap(1, 1),
                                   tmpPlaneRotationMatrixMap(1, 2),
                                   tmpPlaneRotationMatrixMap(2, 0),
                                   tmpPlaneRotationMatrixMap(2, 1),
                                   tmpPlaneRotationMatrixMap(2, 2));
  
   planeUnitVector = tmpPlaneRotationMatrixMap * planeUnitVector;
  
   double denominator = ROOT::Math::Dot(lineUnitVector, planeUnitVector);
   if (denominator == 0) {
     edm::LogError("RPixPlaneCombinatoryTracking")
         << "Error in RPixPlaneCombinatoryTracking::calculatePointOnDetector -> "
         << "Fitted line and plane are parallel. Removing this track";
     return false;
   }
  
   double distanceFromLinePoint = ROOT::Math::Dot((pointOnPlane - pointOnLine), planeUnitVector) / denominator;
  
   math::Vector<3>::type tmpPlaneLineIntercept = distanceFromLinePoint * lineUnitVector + pointOnLine;
   planeLineIntercept = GlobalPoint(tmpPlaneLineIntercept[0], tmpPlaneLineIntercept[1], tmpPlaneLineIntercept[2]);
  
   return true;
 }

References HLTTauDQMOffline_cfi::denominator, RPixDetTrackFinder::geometry_, CTPPSGeometry::localToGlobal(), DetGeomDesc::rotation(), CTPPSGeometry::sensor(), HLT_2018_cff::track, PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), PV3DBase< T, PVType, FrameType >::z(), and HLTMuonOfflineAnalyzer_cfi::z0.

Referenced by findTracks(), and fitTrack().

◆ factorial()

uint32_t RPixPlaneCombinatoryTracking::factorial ( uint32_t x ) const

inlineprivate

Definition at line 61 of file RPixPlaneCombinatoryTracking.h.

                                               {
     if (x == 0)
       return 1;
     return (x == 1 ? x : x * factorial(x - 1));
   }

References x.

Referenced by initialize().

◆ findTracks()

void RPixPlaneCombinatoryTracking::findTracks ( int run )

overridevirtual

Implements RPixDetTrackFinder.

Definition at line 160 of file RPixPlaneCombinatoryTracking.cc.

                                                      {
   //The loop search for all the possible tracks starting from the one with the smallest chiSquare/NDF
   //The loop stops when the number of planes with recorded hits is less than the minimum number of planes required
   //or if the track with minimum chiSquare found has a chiSquare higher than the maximum required
  
   while (hitMap_->size() >= trackMinNumberOfPoints_) {
     if (verbosity_ >= 1)
       edm::LogInfo("RPixPlaneCombinatoryTracking") << "Number of plane with hits " << hitMap_->size();
     if (verbosity_ >= 2)
       for (const auto &plane : *hitMap_)
         edm::LogInfo("RPixPlaneCombinatoryTracking")
             << "\tarm " << plane.first.arm() << " station " << plane.first.station() << " rp " << plane.first.rp()
             << " plane " << plane.first.plane() << " : " << plane.second.size();
  
     //I create the map of all the possible combinations of a group of trackMinNumberOfPoints_ points
     //and the key keeps the reference of which planes and which hit numbers form the combination
     PointAndReferenceMap mapOfAllMinRequiredPoint;
     //I produce the map for all cominations of all hits with all trackMinNumberOfPoints_ plane combinations
     mapOfAllMinRequiredPoint = produceAllHitCombination(possiblePlaneCombinations_);
  
     //Fit all the possible combinations with minimum number of planes required and find the track with minimum chi2
     double theMinChiSquaredOverNDF = maximumChi2OverNDF_ + 1.;  //in order to break the loop in case no track is found;
     HitReferences pointMapWithMinChiSquared;
     PointInPlaneList pointsWithMinChiSquared;
     CTPPSPixelLocalTrack bestTrack;
  
     if (verbosity_ >= 2)
       edm::LogInfo("RPixPlaneCombinatoryTracking")
           << "Number of combinations of trackMinNumberOfPoints_ planes " << mapOfAllMinRequiredPoint.size();
     for (const auto &pointsAndRef : mapOfAllMinRequiredPoint) {
       CTPPSPixelLocalTrack tmpTrack = fitTrack(pointsAndRef.second);
       double tmpChiSquaredOverNDF = tmpTrack.chiSquaredOverNDF();
       if (verbosity_ >= 2)
         edm::LogInfo("RPixPlaneCombinatoryTracking") << "ChiSquare of the present track " << tmpChiSquaredOverNDF;
       if (!tmpTrack.isValid() || tmpChiSquaredOverNDF > maximumChi2OverNDF_ || tmpChiSquaredOverNDF == 0.)
         continue;  //validity check
       if (tmpChiSquaredOverNDF < theMinChiSquaredOverNDF) {
         theMinChiSquaredOverNDF = tmpChiSquaredOverNDF;
         pointMapWithMinChiSquared = pointsAndRef.first;
         pointsWithMinChiSquared = pointsAndRef.second;
         bestTrack = tmpTrack;
       }
     }
  
     //The loop on the fit of all tracks is done, the track with minimum chiSquared is found
     // and it is verified that it complies with the maximumChi2OverNDF_ requirement
     if (theMinChiSquaredOverNDF > maximumChi2OverNDF_)
       break;
  
     //The list of planes not included in the minimum chiSquare track is produced.
     std::vector<uint32_t> listOfExcludedPlanes;
     for (const auto &plane : listOfAllPlanes_) {
       CTPPSPixelDetId tmpRpId = romanPotId_;  //in order to avoid to modify the data member
       tmpRpId.setPlane(plane);
       CTPPSPixelDetId planeDetId = tmpRpId;
       if (pointMapWithMinChiSquared.find(planeDetId) == pointMapWithMinChiSquared.end())
         listOfExcludedPlanes.push_back(plane);
     }
  
     //I produce all the possible combinations of planes to be added to the track,
     //excluding the case of no other plane added since it has been already fitted.
     PlaneCombinations planePointedHitListCombination;
     for (uint32_t i = 1; i <= listOfExcludedPlanes.size(); ++i) {
       PlaneCombinations tmpPlaneCombination;
       getPlaneCombinations(listOfExcludedPlanes, i, tmpPlaneCombination);
       for (const auto &combination : tmpPlaneCombination)
         planePointedHitListCombination.push_back(combination);
     }
  
     //I produce all the possible combinations of points to be added to the track
     PointAndReferenceMap mapOfAllPointWithAtLeastBestFitSelected;
     PointAndReferenceMap mapOfPointCombinationToBeAdded;
     mapOfPointCombinationToBeAdded = produceAllHitCombination(planePointedHitListCombination);
     //The found hit combination is added to the hits selected by the best fit;
     for (const auto &element : mapOfPointCombinationToBeAdded) {
       HitReferences newPointMap = pointMapWithMinChiSquared;
       PointInPlaneList newPoints = pointsWithMinChiSquared;
       for (const auto &pointRef : element.first)
         newPointMap[pointRef.first] = pointRef.second;  //add the new point reference
       for (const auto &point : element.second)
         newPoints.push_back(point);
       mapOfAllPointWithAtLeastBestFitSelected[newPointMap] = newPoints;
     }
  
     //I fit all the possible combination of the minimum plane best fit hits plus hits from the other planes
     if (verbosity_ >= 1)
       edm::LogInfo("RPixPlaneCombinatoryTracking")
           << "Minimum chiSquare over NDF for all the tracks " << theMinChiSquaredOverNDF;
  
     // I look for the tracks with maximum number of points with a chiSquare over NDF smaller than maximumChi2OverNDF_
     // If more than one track fulfill the chiSquare requirement with the same number of points I choose the one with smaller chiSquare
     std::vector<PointAndReferencePair> orderedVectorOfAllPointWithAtLeastBestFitSelected =
         orderCombinationsPerNumberOrPoints(mapOfAllPointWithAtLeastBestFitSelected);
     int currentNumberOfPlanes = 0;
     theMinChiSquaredOverNDF = maximumChi2OverNDF_ + 1.;  //in order to break the loop in case no track is found;
     bool foundTrackWithCurrentNumberOfPlanes = false;
     for (const auto &pointsAndRef : orderedVectorOfAllPointWithAtLeastBestFitSelected) {
       int tmpNumberOfPlanes = pointsAndRef.second.size();
       // If a valid track has been already found with an higher number of planes the loop stops.
       if (foundTrackWithCurrentNumberOfPlanes && tmpNumberOfPlanes < currentNumberOfPlanes)
         break;
       CTPPSPixelLocalTrack tmpTrack = fitTrack(pointsAndRef.second);
       double tmpChiSquaredOverNDF = tmpTrack.chiSquaredOverNDF();
       if (!tmpTrack.isValid() || tmpChiSquaredOverNDF > maximumChi2OverNDF_ || tmpChiSquaredOverNDF == 0.)
         continue;  //validity check
       if (tmpChiSquaredOverNDF < theMinChiSquaredOverNDF) {
         theMinChiSquaredOverNDF = tmpChiSquaredOverNDF;
         pointMapWithMinChiSquared = pointsAndRef.first;
         bestTrack = tmpTrack;
         currentNumberOfPlanes = tmpNumberOfPlanes;
         foundTrackWithCurrentNumberOfPlanes = true;
       }
     }
  
     if (verbosity_ >= 1)
       edm::LogInfo("RPixPlaneCombinatoryTracking") << "The best track has " << bestTrack.ndf() / 2 + 2;
  
     std::vector<uint32_t> listOfPlaneNotUsedForFit = listOfAllPlanes_;
     //remove the hits belonging to the tracks from the full list of hits
     for (const auto &hitToErase : pointMapWithMinChiSquared) {
       std::map<CTPPSPixelDetId, PointInPlaneList>::iterator hitMapElement = hitMap_->find(hitToErase.first);
       if (hitMapElement == hitMap_->end()) {
         throw cms::Exception("RPixPlaneCombinatoryTracking")
             << "The found tracks has hit belonging to a plane which does not have hits";
       }
       std::vector<uint32_t>::iterator planeIt;
       planeIt = std::find(listOfPlaneNotUsedForFit.begin(), listOfPlaneNotUsedForFit.end(), hitToErase.first.plane());
       listOfPlaneNotUsedForFit.erase(planeIt);
       hitMapElement->second.erase(hitMapElement->second.begin() + hitToErase.second);
       //if the plane at which the hit was erased is empty it is removed from the hit map
       if (hitMapElement->second.empty())
         hitMap_->erase(hitMapElement);
     }
  
     //search for hit on the other planes which may belong to the same track
     //even if they did not contributed to the track
     //in case of multiple hit, the closest one to the track will be considered
     //If a hit is found these will not be erased from the list of all hits
     //If no hit is found, the point on the plane intersecting the track will be saved by a CTPPSPixelFittedRecHit
     //with the isRealHit_ flag set to false
     for (const auto &plane : listOfPlaneNotUsedForFit) {
       CTPPSPixelDetId tmpPlaneId = romanPotId_;  //in order to avoid to modify the data member
       tmpPlaneId.setPlane(plane);
       std::unique_ptr<CTPPSPixelFittedRecHit> fittedRecHit(new CTPPSPixelFittedRecHit());
       GlobalPoint pointOnDet;
       calculatePointOnDetector(&bestTrack, tmpPlaneId, pointOnDet);
  
       if (hitMap_->find(tmpPlaneId) != hitMap_->end()) {
         //I convert the hit search window defined in local coordinated into global
         //This avoids to convert the global plane-line intersection in order not to call the the geometry
         math::Vector<3>::type maxGlobalPointDistance(
             maximumXLocalDistanceFromTrack_, maximumYLocalDistanceFromTrack_, 0.);
  
         DetGeomDesc::RotationMatrix theRotationMatrix = geometry_->sensor(tmpPlaneId)->rotation();
         AlgebraicMatrix33 tmpPlaneRotationMatrixMap;
         theRotationMatrix.GetComponents(tmpPlaneRotationMatrixMap(0, 0),
                                         tmpPlaneRotationMatrixMap(0, 1),
                                         tmpPlaneRotationMatrixMap(0, 2),
                                         tmpPlaneRotationMatrixMap(1, 0),
                                         tmpPlaneRotationMatrixMap(1, 1),
                                         tmpPlaneRotationMatrixMap(1, 2),
                                         tmpPlaneRotationMatrixMap(2, 0),
                                         tmpPlaneRotationMatrixMap(2, 1),
                                         tmpPlaneRotationMatrixMap(2, 2));
  
         maxGlobalPointDistance = tmpPlaneRotationMatrixMap * maxGlobalPointDistance;
         //I avoid the Sqrt since it will not be saved
         double maximumXdistance = maxGlobalPointDistance[0] * maxGlobalPointDistance[0];
         double maximumYdistance = maxGlobalPointDistance[1] * maxGlobalPointDistance[1];
         // to be sure that the first min distance is from a real point
         double minimumDistance = 1. + maximumXdistance + maximumYdistance;
         for (const auto &hit : (*hitMap_)[tmpPlaneId]) {
           double xResidual = hit.globalPoint.x() - pointOnDet.x();
           double yResidual = hit.globalPoint.y() - pointOnDet.y();
           double xDistance = xResidual * xResidual;
           double yDistance = yResidual * yResidual;
           double distance = xDistance + yDistance;
           if (xDistance < maximumXdistance && yDistance < maximumYdistance && distance < minimumDistance) {
             LocalPoint residuals(xResidual, yResidual, 0.);
             math::Error<3>::type globalError = hit.globalError;
             LocalPoint pulls(xResidual / std::sqrt(globalError[0][0]), yResidual / std::sqrt(globalError[1][1]), 0.);
             fittedRecHit.reset(new CTPPSPixelFittedRecHit(hit.recHit, pointOnDet, residuals, pulls));
             fittedRecHit->setIsRealHit(true);
           }
         }
       } else {
         LocalPoint fakePoint;
         LocalError fakeError;
         CTPPSPixelRecHit fakeRecHit(fakePoint, fakeError);
         fittedRecHit.reset(new CTPPSPixelFittedRecHit(fakeRecHit, pointOnDet, fakePoint, fakePoint));
       }
  
       bestTrack.addHit(tmpPlaneId, *fittedRecHit);
     }
  
     localTrackVector_.push_back(bestTrack);
  
     int pointForTracking = 0;
     int pointOnTrack = 0;
  
     if (verbosity_ >= 1) {
       for (const auto &planeHits : bestTrack.hits()) {
         for (const auto &fittedhit : planeHits) {
           if (fittedhit.isUsedForFit())
             ++pointForTracking;
           if (fittedhit.isRealHit())
             ++pointOnTrack;
         }
       }
       edm::LogInfo("RPixPlaneCombinatoryTracking")
           << "Best track has " << pointForTracking << " points used for the fit and " << pointOnTrack
           << " points belonging to the track\n";
     }
  
   }  //close of the while loop on all the hits
  
   // recoInfo_ calculation
   // Hardcoded shift periods:
   // Before run 300802: No shift
   // Starting from run 300802: Sec45 St2 Rp3 Pl 0,2,3 ROC 0 shifted.
   // Starting from run 303338: No shift.
   // Starting from run 305169: Sec45 St2 Rp3 Pl 1,3,5 ROC 0 shifted.
   // Starting from run 305965: Sec45 St2 Rp3 Pl 1,3,5 ROC 0 shifted & Sec56 St2 Rp3 Pl 2,4,5 ROC 5 shifted.
   // Starting from run 307083: No shift
  
   // These variables hold the information of the runs when the detector was taking data in 3+3 Mode and which planes were bx-shifted
   // These values will never be changed and the 3+3 Mode will never be used again in the future
   const CTPPSPixelDetId rpId_arm0_st2 = CTPPSPixelDetId(0, 2, 3);
   const CTPPSPixelDetId rpId_arm1_st2 = CTPPSPixelDetId(1, 2, 3);
   static const std::map<unsigned int, std::map<CTPPSPixelDetId, std::vector<bool> > > isPlaneShifted = {
       {0,
        {
            {rpId_arm0_st2, {false, false, false, false, false, false}},  // Shift Period 0 Sec45
            {rpId_arm1_st2, {false, false, false, false, false, false}}   // Shift Period 1 Sec56
        }},
       {300802,
        {
            {rpId_arm0_st2, {true, false, true, true, false, false}},    // Shift Period 1 Sec45
            {rpId_arm1_st2, {false, false, false, false, false, false}}  // Shift Period 1 Sec56
        }},
       {303338,
        {
            {rpId_arm0_st2, {false, false, false, false, false, false}},  // Shift Period 2 Sec45
            {rpId_arm1_st2, {false, false, false, false, false, false}}   // Shift Period 2 Sec56
        }},
       {305169,
        {
            {rpId_arm0_st2, {false, true, false, true, false, true}},    // Shift Period 3 Sec45
            {rpId_arm1_st2, {false, false, false, false, false, false}}  // Shift Period 3 Sec56
        }},
       {305965,
        {
            {rpId_arm0_st2, {false, true, false, true, false, true}},  // Shift Period 4 Sec45
            {rpId_arm1_st2, {false, false, true, false, true, true}}   // Shift Period 4 Sec56
        }},
       {307083,
        {
            {rpId_arm0_st2, {false, false, false, false, false, false}},  // Shift Period 0 Sec45
            {rpId_arm1_st2, {false, false, false, false, false, false}}   // Shift Period 1 Sec56
        }}};                                                              // map< shiftedPeriod, map<DetID,shiftScheme> >
   const auto &shiftStatusInitialRun = std::prev(isPlaneShifted.upper_bound(run));
   unsigned short shiftedROC = 10;
   CTPPSPixelIndices pixelIndices;
  
   // Selecting the shifted ROC
   if (romanPotId_.arm() == 0)
     shiftedROC = 0;
   if (romanPotId_.arm() == 1)
     shiftedROC = 5;
  
   // Loop over found tracks to set recoInfo_
   for (auto &track : localTrackVector_) {
     if (romanPotId_ != rpId_arm0_st2 && romanPotId_ != rpId_arm1_st2) {
       track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::notShiftedRun);
       if (verbosity_ >= 2)
         edm::LogInfo("RPixPlaneCombinatoryTracking") << "Analyzing run: " << run << "\nTrack belongs to Arm "
                                                      << romanPotId_.arm() << " Station " << romanPotId_.station();
  
       continue;
     }
     unsigned short bxShiftedPlanesUsed = 0;
     unsigned short bxNonShiftedPlanesUsed = 0;
     unsigned short hitInShiftedROC = 0;
  
     auto const &fittedHits = track.hits();
     auto const &planeFlags = (shiftStatusInitialRun->second).at(romanPotId_);
  
     for (const auto &planeHits : fittedHits) {
       unsigned short plane = CTPPSPixelDetId(planeHits.detId()).plane();
       for (const auto &hit : planeHits) {
         if (hit.isUsedForFit()) {
           if (pixelIndices.getROCId(hit.minPixelCol(), hit.minPixelRow()) == shiftedROC)
             hitInShiftedROC++;  // Count how many hits are in the shifted ROC
           if (planeFlags.at(plane))
             bxShiftedPlanesUsed++;  // Count how many bx-shifted planes are used
           else if (planeFlags != std::vector<bool>(6, false))
             bxNonShiftedPlanesUsed++;  // Count how many non-bx-shifted planes are used, only if there are shifted planes
         }
       }
     }
  
     // Set recoInfo_ value
     track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
                           invalid);  // Initially setting it as invalid. It has to match one of the following options.
     if (hitInShiftedROC < 3)
       track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::notShiftedRun);
     else {
       if (bxShiftedPlanesUsed == 0 && bxNonShiftedPlanesUsed == 0)
         track.setRecoInfo(
             CTPPSpixelLocalTrackReconstructionInfo::notShiftedRun);  // Default value for runs without bx-shift
       if (bxShiftedPlanesUsed == 3 && bxNonShiftedPlanesUsed == 0)
         track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
                               allShiftedPlanes);  // Track reconstructed in a shifted ROC, only with bx-shifted planes
       if (bxShiftedPlanesUsed == 0 && bxNonShiftedPlanesUsed == 3)
         // Track reconstructed in a shifted ROC, only with non-bx-shifted planes
         track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::noShiftedPlanes);
       if (bxShiftedPlanesUsed > 0 && bxNonShiftedPlanesUsed > 0)
         track.setRecoInfo(CTPPSpixelLocalTrackReconstructionInfo::
                               mixedPlanes);  // Track reconstructed in a shifted ROC, with mixed planes
     }
     if (bxShiftedPlanesUsed + bxNonShiftedPlanesUsed > 6) {
       throw cms::Exception("RPixPlaneCombinatoryTracking") << "Error in RPixPlaneCombinatoryTracking::findTracks -> "
                                                            << "More than six points found for a track, skipping.";
       continue;
     }
     if (track.recoInfo() == CTPPSpixelLocalTrackReconstructionInfo::invalid) {
       throw cms::Exception("RPixPlaneCombinatoryTracking") << "Error in RPixPlaneCombinatoryTracking::findTracks -> "
                                                            << "recoInfo has not been set properly.";
     }
  
     if (verbosity_ >= 2) {
       edm::LogInfo("RPixPlaneCombinatoryTracking")
           << " Track belongs to Arm " << romanPotId_.arm() << " Station " << romanPotId_.station()
           << "\nFirst run with this bx-shift configuration: " << shiftStatusInitialRun->first
           << "\nTrack reconstructed with: " << bxShiftedPlanesUsed << " bx-shifted planes, " << bxNonShiftedPlanesUsed
           << " non-bx-shifted planes, " << hitInShiftedROC << " hits in the bx-shifted ROC"
           << "\nrecoInfo = " << (unsigned short)track.recoInfo();
       if (planeFlags != std::vector<bool>(6, false))
         edm::LogInfo("RPixPlaneCombinatoryTracking") << "The shifted ROC is ROC" << shiftedROC;
     }
   }
   return;
 }

References CTPPSPixelLocalTrack::addHit(), allShiftedPlanes, CTPPSDetId::arm(), calculatePointOnDetector(), CTPPSPixelLocalTrack::chiSquaredOverNDF(), HLT_2018_cff::distance, Exception, spr::find(), fitTrack(), RPixDetTrackFinder::geometry_, getPlaneCombinations(), CTPPSPixelIndices::getROCId(), RPixDetTrackFinder::hitMap_, CTPPSPixelLocalTrack::hits(), mps_fire::i, invalid, CTPPSPixelLocalTrack::isValid(), RPixDetTrackFinder::listOfAllPlanes_, RPixDetTrackFinder::localTrackVector_, maximumChi2OverNDF_, maximumXLocalDistanceFromTrack_, maximumYLocalDistanceFromTrack_, mixedPlanes, CTPPSPixelLocalTrack::ndf(), noShiftedPlanes, notShiftedRun, orderCombinationsPerNumberOrPoints(), point, possiblePlaneCombinations_, produceAllHitCombination(), RPixDetTrackFinder::romanPotId_, DetGeomDesc::rotation(), writedatasetfile::run, CTPPSGeometry::sensor(), CTPPSPixelDetId::setPlane(), mathSSE::sqrt(), CTPPSDetId::station(), HLT_2018_cff::track, trackMinNumberOfPoints_, verbosity_, PV3DBase< T, PVType, FrameType >::x(), hit::x, PV3DBase< T, PVType, FrameType >::y(), and hit::y.

◆ fitTrack()

CTPPSPixelLocalTrack RPixPlaneCombinatoryTracking::fitTrack ( PointInPlaneList pointList )

private

Definition at line 506 of file RPixPlaneCombinatoryTracking.cc.

                                                                                       {
   uint32_t const numberOfPlanes = 6;
   math::Vector<2 * numberOfPlanes>::type xyCoordinates;
   math::Error<2 * numberOfPlanes>::type varianceMatrix;
   math::Matrix<2 * numberOfPlanes, 4>::type zMatrix;
  
   //The matrices and vector xyCoordinates, varianceMatrix and varianceMatrix are built from the points
   for (uint32_t iHit = 0; iHit < numberOfPlanes; iHit++) {
     if (iHit < pointList.size()) {
       const auto &globalPoint = pointList[iHit].globalPoint;
       xyCoordinates[2 * iHit] = globalPoint.x();
       xyCoordinates[2 * iHit + 1] = globalPoint.y();
       zMatrix(2 * iHit, 0) = 1.;
       zMatrix(2 * iHit, 2) = globalPoint.z() - z0_;
       zMatrix(2 * iHit + 1, 1) = 1.;
       zMatrix(2 * iHit + 1, 3) = globalPoint.z() - z0_;
  
       AlgebraicMatrix33 globalError = pointList[iHit].globalError;
       varianceMatrix(2 * iHit, 2 * iHit) = globalError(0, 0);
       varianceMatrix(2 * iHit, 2 * iHit + 1) = globalError(0, 1);
       varianceMatrix(2 * iHit + 1, 2 * iHit) = globalError(1, 0);
       varianceMatrix(2 * iHit + 1, 2 * iHit + 1) = globalError(1, 1);
     } else {
       varianceMatrix(2 * iHit, 2 * iHit) = 1.;
       varianceMatrix(2 * iHit + 1, 2 * iHit + 1) = 1.;
     }
   }
  
   //Get real point variance matrix
   if (!varianceMatrix.Invert()) {
     edm::LogError("RPixPlaneCombinatoryTracking") << "Error in RPixPlaneCombinatoryTracking::fitTrack -> "
                                                   << "Point variance matrix is singular, skipping.";
     CTPPSPixelLocalTrack badTrack;
     badTrack.setValid(false);
     return badTrack;
   }
  
   math::Error<4>::type covarianceMatrix = ROOT::Math::SimilarityT(zMatrix, varianceMatrix);
  
   //To have the real parameter covariance matrix, covarianceMatrix needs to be inverted
   if (!covarianceMatrix.Invert()) {
     edm::LogError("RPixPlaneCombinatoryTracking") << "Error in RPixPlaneCombinatoryTracking::fitTrack -> "
                                                   << "Parameter covariance matrix is singular, skipping.";
     CTPPSPixelLocalTrack badTrack;
     badTrack.setValid(false);
     return badTrack;
   }
  
   // track parameters: (x0, y0, tx, ty); x = x0 + tx*(z-z0)
   math::Vector<4>::type zMatrixTransposeTimesVarianceMatrixTimesXyCoordinates =
       ROOT::Math::Transpose(zMatrix) * varianceMatrix * xyCoordinates;
   math::Vector<4>::type parameterVector = covarianceMatrix * zMatrixTransposeTimesVarianceMatrixTimesXyCoordinates;
   math::Vector<2 *numberOfPlanes>::type xyCoordinatesMinusZmatrixTimesParameters =
       xyCoordinates - (zMatrix * parameterVector);
  
   double chiSquare = ROOT::Math::Dot(xyCoordinatesMinusZmatrixTimesParameters,
                                      (varianceMatrix * xyCoordinatesMinusZmatrixTimesParameters));
  
   CTPPSPixelLocalTrack goodTrack(z0_, parameterVector, covarianceMatrix, chiSquare);
   goodTrack.setValid(true);
  
   for (const auto &hit : pointList) {
     const auto &globalPoint = hit.globalPoint;
     GlobalPoint pointOnDet;
     bool foundPoint = calculatePointOnDetector(&goodTrack, hit.detId, pointOnDet);
     if (!foundPoint) {
       CTPPSPixelLocalTrack badTrack;
       badTrack.setValid(false);
       return badTrack;
     }
     double xResidual = globalPoint.x() - pointOnDet.x();
     double yResidual = globalPoint.y() - pointOnDet.y();
     LocalPoint residuals(xResidual, yResidual);
  
     math::Error<3>::type globalError(hit.globalError);
     LocalPoint pulls(xResidual / std::sqrt(globalError(0, 0)), yResidual / std::sqrt(globalError(0, 0)));
  
     CTPPSPixelFittedRecHit fittedRecHit(hit.recHit, pointOnDet, residuals, pulls);
     fittedRecHit.setIsUsedForFit(true);
     goodTrack.addHit(hit.detId, fittedRecHit);
   }
  
   return goodTrack;
 }

References calculatePointOnDetector(), simKBmtfDigis_cfi::chiSquare, spr::goodTrack(), CTPPSPixelFittedRecHit::setIsUsedForFit(), CTPPSPixelLocalTrack::setValid(), mathSSE::sqrt(), PV3DBase< T, PVType, FrameType >::x(), PV3DBase< T, PVType, FrameType >::y(), and RPixDetTrackFinder::z0_.

Referenced by findTracks().

◆ functionForPlaneOrdering()

static bool RPixPlaneCombinatoryTracking::functionForPlaneOrdering	(	PointAndReferencePair	a,
		PointAndReferencePair	b
	)

inlinestaticprivate

Definition at line 56 of file RPixPlaneCombinatoryTracking.h.

                                                                                          {
     return (a.second.size() > b.second.size());
   }

References a, and b.

Referenced by orderCombinationsPerNumberOrPoints().

◆ getHitCombinations()

void RPixPlaneCombinatoryTracking::getHitCombinations	(	const std::map< CTPPSPixelDetId, PointInPlaneList > &	mapOfAllHits,
		std::map< CTPPSPixelDetId, PointInPlaneList >::iterator	mapIterator,
		HitReferences	tmpHitPlaneMap,
		const PointInPlaneList &	tmpHitVector,
		PointAndReferenceMap &	outputMap
	)

private

Definition at line 87 of file RPixPlaneCombinatoryTracking.cc.

                                                                                        {
   //At this point I selected one hit per plane
   if (mapIterator == mapOfAllHits.end()) {
     outputMap[tmpHitPlaneMap] = tmpHitVector;
     return;
   }
   for (size_t i = 0; i < mapIterator->second.size(); i++) {
     HitReferences newHitPlaneMap = tmpHitPlaneMap;
     newHitPlaneMap[mapIterator->first] = i;
     PointInPlaneList newVector = tmpHitVector;
     newVector.push_back(mapIterator->second[i]);
     std::map<CTPPSPixelDetId, PointInPlaneList>::iterator tmpMapIterator = mapIterator;
     getHitCombinations(mapOfAllHits, ++tmpMapIterator, newHitPlaneMap, newVector, outputMap);
   }
 }

References mps_fire::i.

Referenced by produceAllHitCombination().

◆ getPlaneCombinations()

void RPixPlaneCombinatoryTracking::getPlaneCombinations	(	const std::vector< uint32_t > &	inputPlaneList,
		uint32_t	numberToExtract,
		PlaneCombinations &	planeCombinations
	)		const

private

Definition at line 61 of file RPixPlaneCombinatoryTracking.cc.

                                                                                                     {
   uint32_t numberOfPlanes = inputPlaneList.size();
   std::string bitmask(numberToExtract, 1);  // numberToExtract leading 1's
   bitmask.resize(numberOfPlanes, 0);        // numberOfPlanes-numberToExtract trailing 0's
   planeCombinations.clear();
  
   // store the combination and permute bitmask
   do {
     planeCombinations.push_back(std::vector<uint32_t>());
     for (uint32_t i = 0; i < numberOfPlanes; ++i) {  // [0..numberOfPlanes-1] integers
       if (bitmask[i])
         planeCombinations.back().push_back(inputPlaneList.at(i));
     }
   } while (std::prev_permutation(bitmask.begin(), bitmask.end()));
  
   return;
 }

References mps_fire::i, and AlCaHLTBitMon_QueryRunRegistry::string.

Referenced by findTracks(), and initialize().

◆ initialize()

void RPixPlaneCombinatoryTracking::initialize ( void )

overridevirtual

Implements RPixDetTrackFinder.

Definition at line 38 of file RPixPlaneCombinatoryTracking.cc.

                                               {
   uint32_t numberOfCombinations =
       factorial(numberOfPlanesPerPot_) /
       (factorial(numberOfPlanesPerPot_ - trackMinNumberOfPoints_) * factorial(trackMinNumberOfPoints_));
   if (verbosity_ >= 2)
     edm::LogInfo("RPixPlaneCombinatoryTracking") << "Number of combinations = " << numberOfCombinations;
   possiblePlaneCombinations_.reserve(numberOfCombinations);
  
   getPlaneCombinations(listOfAllPlanes_, trackMinNumberOfPoints_, possiblePlaneCombinations_);
  
   if (verbosity_ >= 2) {
     for (const auto &vec : possiblePlaneCombinations_) {
       for (const auto &num : vec) {
         edm::LogInfo("RPixPlaneCombinatoryTracking") << num << " - ";
       }
       edm::LogInfo("RPixPlaneCombinatoryTracking");
     }
   }
 }

References factorial(), getPlaneCombinations(), RPixDetTrackFinder::listOfAllPlanes_, EgammaValidation_cff::num, RPixDetTrackFinder::numberOfPlanesPerPot_, possiblePlaneCombinations_, trackMinNumberOfPoints_, and verbosity_.

◆ orderCombinationsPerNumberOrPoints()

std::vector< RPixPlaneCombinatoryTracking::PointAndReferencePair > RPixPlaneCombinatoryTracking::orderCombinationsPerNumberOrPoints ( PointAndReferenceMap inputMap )

private

Definition at line 643 of file RPixPlaneCombinatoryTracking.cc.

                                                                                               {
   std::vector<PointAndReferencePair> sortedVector(inputMap.begin(), inputMap.end());
   std::sort(sortedVector.begin(), sortedVector.end(), functionForPlaneOrdering);
  
   return sortedVector;
 }

References functionForPlaneOrdering().

Referenced by findTracks().

◆ produceAllHitCombination()

RPixPlaneCombinatoryTracking::PointAndReferenceMap RPixPlaneCombinatoryTracking::produceAllHitCombination ( PlaneCombinations inputPlaneCombination )

private

Definition at line 109 of file RPixPlaneCombinatoryTracking.cc.

                                              {
   PointAndReferenceMap mapOfAllPoints;
   CTPPSPixelDetId tmpRpId = romanPotId_;  //in order to avoid to modify the data member
  
   if (verbosity_ >= 2)
     edm::LogInfo("RPixPlaneCombinatoryTracking") << "Searching for all combinations...";
   //Loop on all the plane combinations
   for (const auto &planeCombination : inputPlaneCombination) {
     std::map<CTPPSPixelDetId, PointInPlaneList> selectedCombinationHitOnPlane;
     bool allPlaneAsHits = true;
  
     //Loop on all the possible combinations
     //In this loop the selectedCombinationHitOnPlane is filled
     //and cases in which one of the selected plane is empty are skipped
     for (const auto &plane : planeCombination) {
       tmpRpId.setPlane(plane);
       CTPPSPixelDetId planeDetId = tmpRpId;
       if (hitMap_->find(planeDetId) == hitMap_->end()) {
         if (verbosity_ >= 2)
           edm::LogInfo("RPixPlaneCombinatoryTracking")
               << "No data on arm " << planeDetId.arm() << " station " << planeDetId.station() << " rp "
               << planeDetId.rp() << " plane " << planeDetId.plane();
         allPlaneAsHits = false;
         break;
       }
       if (selectedCombinationHitOnPlane.find(planeDetId) != selectedCombinationHitOnPlane.end()) {
         throw cms::Exception("RPixPlaneCombinatoryTracking")
             << "selectedCombinationHitOnPlane contains already detId " << planeDetId
             << "Error in the algorithm which created all the possible plane combinations";
       }
       selectedCombinationHitOnPlane[planeDetId] = (*hitMap_)[planeDetId];
     }
     if (!allPlaneAsHits)
       continue;
  
     //I add the all the hit combinations to the full list of plane combinations
     auto mapIterator = selectedCombinationHitOnPlane.begin();
     HitReferences tmpHitPlaneMap;   //empty map of plane id and hit number needed the getHitCombinations algorithm
     PointInPlaneList tmpHitVector;  //empty vector of hits needed for the getHitCombinations algorithm
     getHitCombinations(selectedCombinationHitOnPlane, mapIterator, tmpHitPlaneMap, tmpHitVector, mapOfAllPoints);
     if (verbosity_ >= 2)
       edm::LogInfo("RPixPlaneCombinatoryTracking") << "Number of possible tracks " << mapOfAllPoints.size();
  
   }  //end of loops on all the combinations
  
   return mapOfAllPoints;
 }