RPixPlaneCombinatoryTracking Class Reference

#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 ()
std::vector< CTPPSPixelLocalTrack > const &	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 17 of file RPixPlaneCombinatoryTracking.cc.

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

    : 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";
  }
}

~RPixPlaneCombinatoryTracking()

RPixPlaneCombinatoryTracking::~RPixPlaneCombinatoryTracking ( )

override

Definition at line 35 of file RPixPlaneCombinatoryTracking.cc.

References possiblePlaneCombinations_.

{ possiblePlaneCombinations_.clear(); }

Member Function Documentation

calculatePointOnDetector()

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

private

Definition at line 647 of file RPixPlaneCombinatoryTracking.cc.

References bTagMiniDQMDeepCSV::denominator, RPixDetTrackFinder::geometry_, CTPPSGeometry::localToGlobal(), DetGeomDesc::rotation(), CTPPSGeometry::sensor(), HLT_2022v15_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().

                                                                                             {
  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;
}

factorial()

uint32_t RPixPlaneCombinatoryTracking::factorial ( uint32_t  x ) const

inlineprivate

Definition at line 61 of file RPixPlaneCombinatoryTracking.h.

References x.

Referenced by initialize().

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

findTracks()

void RPixPlaneCombinatoryTracking::findTracks ( int  run )

overridevirtual

Implements RPixDetTrackFinder.

Definition at line 161 of file RPixPlaneCombinatoryTracking.cc.

References CTPPSPixelLocalTrack::addHit(), allShiftedPlanes, CTPPSDetId::arm(), calculatePointOnDetector(), CTPPSPixelLocalTrack::chiSquaredOverNDF(), HLT_2022v15_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(), CTPPSGeometry::rpTranslation(), writedatasetfile::run, CTPPSGeometry::sensor(), CTPPSPixelFittedRecHit::setIsRealHit(), CTPPSPixelDetId::setPlane(), mathSSE::sqrt(), CTPPSDetId::station(), HLT_2022v15_cff::track, trackMinNumberOfPoints_, findQualityFiles::v, verbosity_, PV3DBase< T, PVType, FrameType >::x(), hit::x, PV3DBase< T, PVType, FrameType >::y(), hit::y, z, hit::z, and HLTMuonOfflineAnalyzer_cfi::z0.

                                                     {
  //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

  // bad Pot patch 45-220-fr 2022 -- beginning
  // check number of hits in road
  unsigned int hitNum = 0;
  for (const auto &plane : *hitMap_) {
    hitNum += plane.second.size();
  }

  if (hitMap_->size() == 2 && hitNum == 2) {  // look for roads with 2 hits in 2 different planes

    GlobalPoint hit[2];
    PointInPlaneList pIPL;

    unsigned int i = 0;
    for (const auto &plane : *hitMap_) {
      if (plane.second.size() > 1)
        break;  // only 1 hit per plane per road allowed
      GlobalPoint gP(plane.second[0].globalPoint.x(), plane.second[0].globalPoint.y(), plane.second[0].globalPoint.z());
      hit[i] = gP;
      i++;
      pIPL.push_back(plane.second[0]);
    }

    // calculate 2 point track parameters (no need of fits)
    double tx = (hit[0].x() - hit[1].x()) / (hit[0].z() - hit[1].z());
    double ty = (hit[0].y() - hit[1].y()) / (hit[0].z() - hit[1].z());
    double xat0 = (hit[1].x() * hit[0].z() - hit[0].x() * hit[1].z()) / (hit[0].z() - hit[1].z());
    double yat0 = (hit[1].y() * hit[0].z() - hit[0].y() * hit[1].z()) / (hit[0].z() - hit[1].z());
    double z0 = geometry_->rpTranslation(romanPotId_).z();
    double xatz0 = xat0 + tx * z0;
    double yatz0 = yat0 + ty * z0;

    math::Vector<4>::type parameterVector{xatz0, yatz0, tx, ty};
    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);
    math::Error<4>::type covarianceMatrix(v);

    CTPPSPixelLocalTrack track(z0, parameterVector, covarianceMatrix, 0);

    // save used points into track
    for (const auto &hhit : pIPL) {
      GlobalPoint pOD(hhit.globalPoint.x(), hhit.globalPoint.y(), hhit.globalPoint.z());
      LocalPoint res(0, 0);
      LocalPoint pulls(0, 0);
      CTPPSPixelFittedRecHit usedRecHit(hhit.recHit, pOD, res, pulls);
      usedRecHit.setIsRealHit(true);
      track.addHit(hhit.detId, usedRecHit);
    }

    // save track in collection
    localTrackVector_.push_back(track);
  }
  // bad Pot patch 45-220-fr 2022 -- end

  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 = std::make_unique<CTPPSPixelFittedRecHit>(hit.recHit, pointOnDet, residuals, pulls);
            fittedRecHit->setIsRealHit(true);
          }
        }
      } else {
        LocalPoint fakePoint;
        LocalError fakeError;
        CTPPSPixelRecHit fakeRecHit(fakePoint, fakeError);
        fittedRecHit = std::make_unique<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;
}

fitTrack()

CTPPSPixelLocalTrack RPixPlaneCombinatoryTracking::fitTrack ( PointInPlaneList  pointList )

private

Definition at line 559 of file RPixPlaneCombinatoryTracking.cc.

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().

                                                                                      {
  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;
}

functionForPlaneOrdering()

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

inlinestaticprivate

Definition at line 56 of file RPixPlaneCombinatoryTracking.h.

References a, and b.

Referenced by orderCombinationsPerNumberOrPoints().

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

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 88 of file RPixPlaneCombinatoryTracking.cc.

References mps_fire::i.

Referenced by produceAllHitCombination().

                                                                                       {
  //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);
  }
}

getPlaneCombinations()

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

private

Definition at line 62 of file RPixPlaneCombinatoryTracking.cc.

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

Referenced by findTracks(), and initialize().

                                                                                                    {
  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.emplace_back();
    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;
}

initialize()

void RPixPlaneCombinatoryTracking::initialize ( void  )

overridevirtual

Implements RPixDetTrackFinder.

Definition at line 39 of file RPixPlaneCombinatoryTracking.cc.

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

                                              {
  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");
    }
  }
}

orderCombinationsPerNumberOrPoints()

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

private

Definition at line 696 of file RPixPlaneCombinatoryTracking.cc.

References functionForPlaneOrdering(), and jetsAK4_CHS_cff::sort.

Referenced by findTracks().

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

  return sortedVector;
}

produceAllHitCombination()

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

private

Definition at line 110 of file RPixPlaneCombinatoryTracking.cc.

References CTPPSDetId::arm(), Exception, getHitCombinations(), RPixDetTrackFinder::hitMap_, CTPPSPixelDetId::plane(), RPixDetTrackFinder::romanPotId_, CTPPSDetId::rp(), CTPPSPixelDetId::setPlane(), CTPPSDetId::station(), and verbosity_.

Referenced by findTracks().

                                             {
  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;
}

Member Data Documentation

maximumChi2OverNDF_

double RPixPlaneCombinatoryTracking::maximumChi2OverNDF_

private

Definition at line 40 of file RPixPlaneCombinatoryTracking.h.

Referenced by findTracks(), and RPixPlaneCombinatoryTracking().

maximumXLocalDistanceFromTrack_

double RPixPlaneCombinatoryTracking::maximumXLocalDistanceFromTrack_

private

Definition at line 41 of file RPixPlaneCombinatoryTracking.h.

Referenced by findTracks(), and RPixPlaneCombinatoryTracking().

maximumYLocalDistanceFromTrack_

double RPixPlaneCombinatoryTracking::maximumYLocalDistanceFromTrack_

private

Definition at line 42 of file RPixPlaneCombinatoryTracking.h.

Referenced by findTracks(), and RPixPlaneCombinatoryTracking().

possiblePlaneCombinations_

PlaneCombinations RPixPlaneCombinatoryTracking::possiblePlaneCombinations_

private

Definition at line 43 of file RPixPlaneCombinatoryTracking.h.

Referenced by findTracks(), initialize(), and ~RPixPlaneCombinatoryTracking().

trackMinNumberOfPoints_

uint32_t RPixPlaneCombinatoryTracking::trackMinNumberOfPoints_

private

Definition at line 39 of file RPixPlaneCombinatoryTracking.h.

Referenced by findTracks(), initialize(), and RPixPlaneCombinatoryTracking().

verbosity_

int RPixPlaneCombinatoryTracking::verbosity_

private

Definition at line 38 of file RPixPlaneCombinatoryTracking.h.

Referenced by findTracks(), initialize(), produceAllHitCombination(), and RPixPlaneCombinatoryTracking().