Merge pull request #38553 from CTPPS/badPotReco2

PPS reco patch for bad pot

RecoPPS/Local/interface/RPixDetPatternFinder.h

@@ -38,7 +38,7 @@ class RPixDetPatternFinder {
   typedef std::vector<PointInPlane> Road;
 
   void setHits(const edm::DetSetVector<CTPPSPixelRecHit> *hitVector) { hitVector_ = hitVector; }
-  virtual void findPattern() = 0;
+  virtual void findPattern(bool isbadpot) = 0;
   void clear() { patternVector_.clear(); }
   std::vector<Road> const &getPatterns() const { return patternVector_; }
   void setGeometry(const CTPPSGeometry *geometry) { geometry_ = geometry; }

RecoPPS/Local/interface/RPixRoadFinder.h

@@ -35,13 +35,15 @@ class RPixRoadFinder : public RPixDetPatternFinder {
 public:
   explicit RPixRoadFinder(const edm::ParameterSet &param);
   ~RPixRoadFinder() override;
-  void findPattern() override;
+  void findPattern(bool isbadpot) override;
 
 private:
   int verbosity_;
   double roadRadius_;
   unsigned int minRoadSize_;
   unsigned int maxRoadSize_;
+  double roadRadiusBadPot_;
+  //  bool isBadPot_;
   void run(const edm::DetSetVector<CTPPSPixelRecHit> &input, const CTPPSGeometry &geometry, std::vector<Road> &roads);
 };
 

RecoPPS/Local/plugins/CTPPSPixelLocalTrackProducer.cc

@@ -38,6 +38,15 @@
 #include "RecoPPS/Local/interface/RPixRoadFinder.h"
 #include "RecoPPS/Local/interface/RPixPlaneCombinatoryTracking.h"
 
+#include "CondFormats/PPSObjects/interface/CTPPSPixelAnalysisMask.h"
+#include "CondFormats/DataRecord/interface/CTPPSPixelAnalysisMaskRcd.h"
+
+namespace {
+  constexpr int rocMask = 0xE000;
+  constexpr int rocOffset = 13;
+  constexpr int rocSizeInPixels = 4160;
+}  // namespace
+
 class CTPPSPixelLocalTrackProducer : public edm::stream::EDProducer<> {
 public:
   explicit CTPPSPixelLocalTrackProducer(const edm::ParameterSet &parameterSet);
@@ -58,13 +67,16 @@ class CTPPSPixelLocalTrackProducer : public edm::stream::EDProducer<> {
   edm::EDGetTokenT<edm::DetSetVector<CTPPSPixelRecHit>> tokenCTPPSPixelRecHit_;
   edm::ESGetToken<CTPPSGeometry, VeryForwardRealGeometryRecord> tokenCTPPSGeometry_;
   edm::ESWatcher<VeryForwardRealGeometryRecord> geometryWatcher_;
+
+  edm::ESGetToken<CTPPSPixelAnalysisMask, CTPPSPixelAnalysisMaskRcd> tokenCTPPSPixelAnalysisMask_;
+
   uint32_t numberOfPlanesPerPot_;
   std::vector<uint32_t> listOfAllPlanes_;
 
   std::unique_ptr<RPixDetPatternFinder> patternFinder_;
   std::unique_ptr<RPixDetTrackFinder> trackFinder_;
 
-  void run(const edm::DetSetVector<CTPPSPixelRecHit> &input, edm::DetSetVector<CTPPSPixelLocalTrack> &output);
+  //  void run(const edm::DetSetVector<CTPPSPixelRecHit> &input, edm::DetSetVector<CTPPSPixelLocalTrack> &output);
 };
 
 //------------------------------------------------------------------------------------------------//
@@ -106,6 +118,7 @@ CTPPSPixelLocalTrackProducer::CTPPSPixelLocalTrackProducer(const edm::ParameterS
 
   tokenCTPPSPixelRecHit_ = consumes<edm::DetSetVector<CTPPSPixelRecHit>>(inputTag_);
   tokenCTPPSGeometry_ = esConsumes<CTPPSGeometry, VeryForwardRealGeometryRecord>();
+  tokenCTPPSPixelAnalysisMask_ = esConsumes<CTPPSPixelAnalysisMask, CTPPSPixelAnalysisMaskRcd>();
 
   produces<edm::DetSetVector<CTPPSPixelLocalTrack>>();
 }
@@ -143,6 +156,9 @@ void CTPPSPixelLocalTrackProducer::fillDescriptions(edm::ConfigurationDescriptio
   desc.add<double>("roadRadius", 1.0)->setComment("radius of pattern search window");
   desc.add<int>("minRoadSize", 3)->setComment("minimum number of points in a pattern");
   desc.add<int>("maxRoadSize", 20)->setComment("maximum number of points in a pattern");
+  //parameters for bad pot reconstruction patch 45-220-fr 2022
+  desc.add<double>("roadRadiusBadPot", 0.5)->setComment("radius of pattern search window for bad Pot");
+  //  desc.add<bool>("isBadPot", true)->setComment("flag to enable road search for bad pot");
 
   descriptions.add("ctppsPixelLocalTracks", desc);
 }
@@ -161,6 +177,32 @@ void CTPPSPixelLocalTrackProducer::produce(edm::Event &iEvent, const edm::EventS
   const CTPPSGeometry &geometry = *geometryHandler;
   geometryWatcher_.check(iSetup);
 
+  // get mask
+  bool isBadPot_45_220 = false;
+  if (!recHits->empty()) {
+    const auto &mask = iSetup.getData(tokenCTPPSPixelAnalysisMask_);
+
+    // Read Mask checking if 45-220-far is masked as bad and needs special treatment
+    std::map<uint32_t, CTPPSPixelROCAnalysisMask> const &maschera = mask.analysisMask;
+
+    bool mask_45_220[6][6] = {{false}};
+    for (auto const &det : maschera) {
+      CTPPSPixelDetId detId(det.first);
+      unsigned int rocNum = (det.first & rocMask) >> rocOffset;
+      if (rocNum > 5 || detId.plane() > 5)
+        throw cms::Exception("InvalidRocOrPlaneNumber") << "roc number from mask: " << rocNum;
+
+      if (detId.arm() == 0 && detId.station() == 2 && detId.rp() == 3) {  // pot 45-220-far
+        if (det.second.maskedPixels.size() == rocSizeInPixels) {          // roc fully masked
+          mask_45_220[detId.plane()][rocNum] = true;
+        }
+      }
+    }
+
+    // search for specific pattern that requires special reconstruction (isBadPot)
+    isBadPot_45_220 = mask_45_220[1][4] && mask_45_220[1][5] && mask_45_220[2][4] && mask_45_220[2][5] &&
+                      mask_45_220[3][4] && mask_45_220[3][5] && mask_45_220[4][4] && mask_45_220[4][5];
+  }
   std::vector<CTPPSPixelDetId> listOfPotWithHighOccupancyPlanes;
   std::map<CTPPSPixelDetId, uint32_t> mapHitPerPot;
 
@@ -206,7 +248,7 @@ void CTPPSPixelLocalTrackProducer::produce(edm::Event &iEvent, const edm::EventS
   patternFinder_->clear();
   patternFinder_->setHits(&recHitVector);
   patternFinder_->setGeometry(&geometry);
-  patternFinder_->findPattern();
+  patternFinder_->findPattern(isBadPot_45_220);
   std::vector<RPixDetPatternFinder::Road> patternVector = patternFinder_->getPatterns();
 
   //loop on all the patterns

RecoPPS/Local/python/ctppsPixelLocalReconstruction_cff.py

@@ -9,6 +9,11 @@
 # local track producer
 from RecoPPS.Local.ctppsPixelLocalTracks_cfi import ctppsPixelLocalTracks
 
+#from Configuration.Eras.Modifier_ctpps_2016_cff import ctpps_2016
+#from Configuration.Eras.Modifier_ctpps_2017_cff import ctpps_2017
+#from Configuration.Eras.Modifier_ctpps_2018_cff import ctpps_2018
+#(ctpps_2016 | ctpps_2017 | ctpps_2018).toModify(ctppsPixelLocalTracks, isBadPot = cms.bool(False))
+
 ctppsPixelLocalReconstructionTask = cms.Task(
     ctppsPixelClusters,ctppsPixelRecHits,ctppsPixelLocalTracks
 )

RecoPPS/Local/src/RPixPlaneCombinatoryTracking.cc

@@ -163,6 +163,58 @@ void RPixPlaneCombinatoryTracking::findTracks(int run) {
   //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();

RecoPPS/Local/src/RPixRoadFinder.cc

@@ -20,6 +20,8 @@ RPixRoadFinder::RPixRoadFinder(edm::ParameterSet const& parameterSet) : RPixDetP
   roadRadius_ = parameterSet.getParameter<double>("roadRadius");
   minRoadSize_ = parameterSet.getParameter<int>("minRoadSize");
   maxRoadSize_ = parameterSet.getParameter<int>("maxRoadSize");
+  roadRadiusBadPot_ = parameterSet.getParameter<double>("roadRadiusBadPot");
+  //  isBadPot_ = parameterSet.getParameter<bool>("isBadPot");
 }
 
 //------------------------------------------------------------------------------------------------//
@@ -28,10 +30,13 @@ RPixRoadFinder::~RPixRoadFinder() {}
 
 //------------------------------------------------------------------------------------------------//
 
-void RPixRoadFinder::findPattern() {
+void RPixRoadFinder::findPattern(bool isBadPot) {
   Road temp_all_hits;
   temp_all_hits.clear();
 
+  Road temp_all_hits_badPot;
+  temp_all_hits_badPot.clear();
+
   // convert local hit sto global and push them to a vector
   for (const auto& ds_rh2 : *hitVector_) {
     const auto myid = CTPPSPixelDetId(ds_rh2.id);
@@ -64,6 +69,14 @@ void RPixRoadFinder::findPattern() {
                                       theRotationTMatrix(2, 2));
 
       math::Error<3>::type globalError = ROOT::Math::SimilarityT(theRotationTMatrix, localError);
+
+      // create new collection for planes 0 and 5 of pot 45-220-fr
+
+      if (isBadPot == true && myid.arm() == 0 && myid.station() == 2 && localV.x() > 0 &&
+          (myid.plane() == 0 || myid.plane() == 5)) {  // 45-220-far
+
+        temp_all_hits_badPot.emplace_back(PointInPlane{globalV, globalError, it_rh, myid});
+      }
       temp_all_hits.emplace_back(PointInPlane{globalV, globalError, it_rh, myid});
     }
   }
@@ -102,4 +115,31 @@ void RPixRoadFinder::findPattern() {
       patternVector_.push_back(temp_road);
   }
   // end of algorithm
+
+  // badPot algorithm
+  Road::iterator it_gh1_bP = temp_all_hits_badPot.begin();
+  Road::iterator it_gh2_bP;
+
+  while (it_gh1_bP != temp_all_hits_badPot.end() && temp_all_hits_badPot.size() >= 2) {
+    Road temp_road;
+
+    it_gh2_bP = it_gh1_bP;
+
+    const auto currPoint = it_gh1_bP->globalPoint;
+
+    while (it_gh2_bP != temp_all_hits_badPot.end()) {
+      const auto subtraction = currPoint - it_gh2_bP->globalPoint;
+
+      if (subtraction.Rho() < roadRadiusBadPot_) {
+        temp_road.push_back(*it_gh2_bP);
+        temp_all_hits_badPot.erase(it_gh2_bP);
+      } else {
+        ++it_gh2_bP;
+      }
+    }
+
+    if (temp_road.size() == 2) {  // look for isolated tracks
+      patternVector_.push_back(temp_road);
+    }
+  }
 }

Validation/CTPPS/test/simu/template_cfg.py

@@ -7,6 +7,16 @@
 import Validation.CTPPS.simu_config.year_$CONFIG_cff as config
 process.load("Validation.CTPPS.simu_config.year_$CONFIG_cff")
 
+process.load("CondCore.CondDB.CondDB_cfi")
+process.CondDB.connect = 'frontier://FrontierProd/CMS_CONDITIONS'
+process.PoolDBESSource = cms.ESSource("PoolDBESSource",
+    process.CondDB,
+    toGet = cms.VPSet(cms.PSet(
+        record = cms.string('CTPPSPixelAnalysisMaskRcd'),
+        tag = cms.string("CTPPSPixelAnalysisMask_Run3_v1_hlt")
+    ))
+)
+
 # minimal logger settings
 process.MessageLogger = cms.Service("MessageLogger",
   statistics = cms.untracked.vstring(),
@@ -63,4 +73,4 @@
   * process.ctppsLHCInfoPlotter
   * process.ctppsTrackDistributionPlotter
   * process.ctppsProtonReconstructionPlotter
-)
+)