Local/src/RPixRoadFinder.cc

0001
0002 #include "RecoPPS/Local/interface/RPixRoadFinder.h"
0003
0004 // user include files
0005 #include "FWCore/Framework/interface/Frameworkfwd.h"
0006
0007 #include "TMath.h"
0008 #include "DataFormats/Math/interface/Error.h"
0009 #include "DataFormats/Math/interface/AlgebraicROOTObjects.h"
0010
0011 #include <vector>
0012 #include <memory>
0013 #include <string>
0014 #include <iostream>
0015
0016 //------------------------------------------------------------------------------------------------//
0017
0018 RPixRoadFinder::RPixRoadFinder(edm::ParameterSet const& parameterSet) : RPixDetPatternFinder(parameterSet) {
0019   verbosity_ = parameterSet.getUntrackedParameter<int>("verbosity");
0020   roadRadius_ = parameterSet.getParameter<double>("roadRadius");
0021   minRoadSize_ = parameterSet.getParameter<int>("minRoadSize");
0022   maxRoadSize_ = parameterSet.getParameter<int>("maxRoadSize");
0023   roadRadiusBadPot_ = parameterSet.getParameter<double>("roadRadiusBadPot");
0024 }
0025
0026 //------------------------------------------------------------------------------------------------//
0027
0028 RPixRoadFinder::~RPixRoadFinder() {}
0029
0030 //------------------------------------------------------------------------------------------------//
0031
0032 void RPixRoadFinder::findPattern(bool* is2PlanePot) {
0033   Road temp_all_hits;
0034   Road temp_all_hits_2PlanePot[4];
0035
0036   // convert local hit sto global and push them to a vector
0037   for (const auto& ds_rh2 : *hitVector_) {
0038     const auto myid = CTPPSPixelDetId(ds_rh2.id);
0039     for (const auto& it_rh : ds_rh2.data) {
0040       CTPPSGeometry::Vector localV(it_rh.point().x(), it_rh.point().y(), it_rh.point().z());
0041       const auto& globalV = geometry_->localToGlobal(ds_rh2.id, localV);
0042       math::Error<3>::type localError;
0043       localError[0][0] = it_rh.error().xx();
0044       localError[0][1] = it_rh.error().xy();
0045       localError[0][2] = 0.;
0046       localError[1][0] = it_rh.error().xy();
0047       localError[1][1] = it_rh.error().yy();
0048       localError[1][2] = 0.;
0049       localError[2][0] = 0.;
0050       localError[2][1] = 0.;
0051       localError[2][2] = 0.;
0052       if (verbosity_ > 2)
0053         edm::LogInfo("RPixRoadFinder") << "Hits = " << ds_rh2.data.size();
0054
0055       DetGeomDesc::RotationMatrix theRotationMatrix = geometry_->sensor(myid)->rotation();
0056       AlgebraicMatrix33 theRotationTMatrix;
0057       theRotationMatrix.GetComponents(theRotationTMatrix(0, 0),
0058                                       theRotationTMatrix(0, 1),
0059                                       theRotationTMatrix(0, 2),
0060                                       theRotationTMatrix(1, 0),
0061                                       theRotationTMatrix(1, 1),
0062                                       theRotationTMatrix(1, 2),
0063                                       theRotationTMatrix(2, 0),
0064                                       theRotationTMatrix(2, 1),
0065                                       theRotationTMatrix(2, 2));
0066
0067       math::Error<3>::type globalError = ROOT::Math::SimilarityT(theRotationTMatrix, localError);
0068
0069       // create new collections for bad (2 planes) and good pots
0070
0071       if (is2PlanePot[0] == true && myid.arm() == 0 && myid.station() == 2) {  // 45-220
0072         temp_all_hits_2PlanePot[0].emplace_back(PointInPlane{globalV, globalError, it_rh, myid});
0073       } else if (is2PlanePot[1] == true && myid.arm() == 0 && myid.station() == 0) {  // 45-210
0074         temp_all_hits_2PlanePot[1].emplace_back(PointInPlane{globalV, globalError, it_rh, myid});
0075       } else if (is2PlanePot[2] == true && myid.arm() == 1 && myid.station() == 0) {  // 56-210
0076         temp_all_hits_2PlanePot[2].emplace_back(PointInPlane{globalV, globalError, it_rh, myid});
0077       } else if (is2PlanePot[3] == true && myid.arm() == 1 && myid.station() == 2) {  // 56-220
0078         temp_all_hits_2PlanePot[3].emplace_back(PointInPlane{globalV, globalError, it_rh, myid});
0079       } else {
0080         temp_all_hits.emplace_back(PointInPlane{globalV, globalError, it_rh, myid});
0081       }
0082     }
0083   }
0084
0085   Road::iterator it_gh1 = temp_all_hits.begin();
0086   Road::iterator it_gh2;
0087
0088   patternVector_.clear();
0089
0090   //look for points near wrt each other
0091   // starting algorithm
0092   while (it_gh1 != temp_all_hits.end() && temp_all_hits.size() >= minRoadSize_) {
0093     Road temp_road;
0094
0095     it_gh2 = it_gh1;
0096
0097     const auto currPoint = it_gh1->globalPoint;
0098     CTPPSPixelDetId currDet = CTPPSPixelDetId(it_gh1->detId);
0099
0100     while (it_gh2 != temp_all_hits.end()) {
0101       bool same_pot = false;
0102       CTPPSPixelDetId tmpGh2Id = CTPPSPixelDetId(it_gh2->detId);
0103       if (currDet.rpId() == tmpGh2Id.rpId())
0104         same_pot = true;
0105       const auto subtraction = currPoint - it_gh2->globalPoint;
0106
0107       if (subtraction.Rho() < roadRadius_ && same_pot) {  /// 1mm
0108         temp_road.push_back(*it_gh2);
0109         temp_all_hits.erase(it_gh2);
0110       } else {
0111         ++it_gh2;
0112       }
0113     }
0114
0115     if (temp_road.size() >= minRoadSize_ && temp_road.size() < maxRoadSize_)
0116       patternVector_.push_back(temp_road);
0117   }
0118   // end of algorithm
0119
0120   // 2PlanePot algorithm
0121
0122   for (unsigned int i = 0; i < 4; i++) {
0123     if (is2PlanePot[i]) {
0124       Road::iterator it_gh1_bP = temp_all_hits_2PlanePot[i].begin();
0125       Road::iterator it_gh2_bP;
0126
0127       while (it_gh1_bP != temp_all_hits_2PlanePot[i].end() && temp_all_hits_2PlanePot[i].size() >= 2) {
0128         Road temp_road;
0129
0130         it_gh2_bP = it_gh1_bP;
0131
0132         const auto currPoint = it_gh1_bP->globalPoint;
0133
0134         while (it_gh2_bP != temp_all_hits_2PlanePot[i].end()) {
0135           const auto subtraction = currPoint - it_gh2_bP->globalPoint;
0136
0137           if (subtraction.Rho() < roadRadiusBadPot_) {
0138             temp_road.push_back(*it_gh2_bP);
0139             temp_all_hits_2PlanePot[i].erase(it_gh2_bP);
0140           } else {
0141             ++it_gh2_bP;
0142           }
0143         }
0144
0145         if (temp_road.size() == 2) {  // look for isolated tracks
0146           patternVector_.push_back(temp_road);
0147         }
0148       }
0149     }
0150   }
0151 }