File indexing completed on 2024-04-06 12:28:35
0001 #include <unordered_map>
0002
0003 #include "CommonTools/Utils/interface/DynArray.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006 #include "FWCore/Framework/interface/Event.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0008 #include "FWCore/Utilities/interface/isFinite.h"
0009 #include "RecoTracker/PixelSeeding/interface/CAHitTripletGenerator.h"
0010 #include "RecoTracker/PixelSeeding/interface/HitTripletGenerator.h"
0011 #include "RecoTracker/PixelSeeding/interface/ThirdHitPredictionFromCircle.h"
0012 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0013
0014 #include "RecoTracker/PixelSeeding/interface/CAGraph.h"
0015 #include "CellularAutomaton.h"
0016
0017 namespace {
0018 template <typename T>
0019 T sqr(T x) {
0020 return x * x;
0021 }
0022 }
0023
0024 using namespace std;
0025
0026 constexpr unsigned int CAHitTripletGenerator::minLayers;
0027
0028 CAHitTripletGenerator::CAHitTripletGenerator(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0029 : theFieldToken(iC.esConsumes()),
0030 extraHitRPhitolerance(cfg.getParameter<double>(
0031 "extraHitRPhitolerance")),
0032 maxChi2(cfg.getParameter<edm::ParameterSet>("maxChi2")),
0033 useBendingCorrection(cfg.getParameter<bool>("useBendingCorrection")),
0034 caThetaCut(cfg.getParameter<double>("CAThetaCut"),
0035 cfg.getParameter<std::vector<edm::ParameterSet>>("CAThetaCut_byTriplets")),
0036 caPhiCut(cfg.getParameter<double>("CAPhiCut"),
0037 cfg.getParameter<std::vector<edm::ParameterSet>>("CAPhiCut_byTriplets")),
0038 caHardPtCut(cfg.getParameter<double>("CAHardPtCut")) {
0039 edm::ParameterSet comparitorPSet = cfg.getParameter<edm::ParameterSet>("SeedComparitorPSet");
0040 std::string comparitorName = comparitorPSet.getParameter<std::string>("ComponentName");
0041 if (comparitorName != "none") {
0042 theComparitor = SeedComparitorFactory::get()->create(comparitorName, comparitorPSet, iC);
0043 }
0044 }
0045
0046 void CAHitTripletGenerator::fillDescriptions(edm::ParameterSetDescription& desc) {
0047 desc.add<double>("extraHitRPhitolerance", 0.06);
0048 desc.add<bool>("useBendingCorrection", false);
0049 desc.add<double>("CAThetaCut", 0.00125);
0050 desc.add<double>("CAPhiCut", 0.1);
0051
0052 edm::ParameterSetDescription validatorCACut;
0053 validatorCACut.add<string>("seedingLayers", "BPix1+BPix2+BPix3");
0054 validatorCACut.add<double>("cut", 0.00125);
0055 std::vector<edm::ParameterSet> defaultCACutVector;
0056 edm::ParameterSet defaultCACut;
0057 defaultCACut.addParameter<string>("seedingLayers", "");
0058 defaultCACut.addParameter<double>("cut", -1.);
0059 defaultCACutVector.push_back(defaultCACut);
0060 desc.addVPSet("CAThetaCut_byTriplets", validatorCACut, defaultCACutVector);
0061 desc.addVPSet("CAPhiCut_byTriplets", validatorCACut, defaultCACutVector);
0062
0063 desc.add<double>("CAHardPtCut", 0);
0064
0065 edm::ParameterSetDescription descMaxChi2;
0066 descMaxChi2.add<double>("pt1", 0.8);
0067 descMaxChi2.add<double>("pt2", 2);
0068 descMaxChi2.add<double>("value1", 50);
0069 descMaxChi2.add<double>("value2", 8);
0070 descMaxChi2.add<bool>("enabled", true);
0071 desc.add<edm::ParameterSetDescription>("maxChi2", descMaxChi2);
0072
0073 edm::ParameterSetDescription descComparitor;
0074 descComparitor.add<std::string>("ComponentName", "none");
0075 descComparitor.setAllowAnything();
0076 desc.add<edm::ParameterSetDescription>("SeedComparitorPSet", descComparitor);
0077 }
0078
0079 void CAHitTripletGenerator::initEvent(const edm::Event& ev, const edm::EventSetup& es) {
0080 if (theComparitor)
0081 theComparitor->init(ev, es);
0082 theField = &es.getData(theFieldToken);
0083 }
0084
0085 namespace {
0086 void createGraphStructure(const SeedingLayerSetsHits& layers, CAGraph& g) {
0087 for (unsigned int i = 0; i < layers.size(); i++) {
0088 for (unsigned int j = 0; j < 3; ++j) {
0089 auto vertexIndex = 0;
0090 auto foundVertex = std::find(g.theLayers.begin(), g.theLayers.end(), layers[i][j].name());
0091 if (foundVertex == g.theLayers.end()) {
0092 g.theLayers.emplace_back(layers[i][j].name(), layers[i][j].detLayer()->seqNum(), layers[i][j].hits().size());
0093 vertexIndex = g.theLayers.size() - 1;
0094 } else {
0095 vertexIndex = foundVertex - g.theLayers.begin();
0096 }
0097 if (j == 0) {
0098 if (std::find(g.theRootLayers.begin(), g.theRootLayers.end(), vertexIndex) == g.theRootLayers.end()) {
0099 g.theRootLayers.emplace_back(vertexIndex);
0100 }
0101 }
0102 }
0103 }
0104 }
0105
0106 void clearGraphStructure(const SeedingLayerSetsHits& layers, CAGraph& g) {
0107 g.theLayerPairs.clear();
0108 for (unsigned int i = 0; i < g.theLayers.size(); i++) {
0109 g.theLayers[i].theInnerLayers.clear();
0110 g.theLayers[i].theInnerLayerPairs.clear();
0111 g.theLayers[i].theOuterLayers.clear();
0112 g.theLayers[i].theOuterLayerPairs.clear();
0113 for (auto& v : g.theLayers[i].isOuterHitOfCell)
0114 v.clear();
0115 }
0116 }
0117
0118 void fillGraph(const SeedingLayerSetsHits& layers,
0119 const IntermediateHitDoublets::RegionLayerSets& regionLayerPairs,
0120 CAGraph& g,
0121 std::vector<const HitDoublets*>& hitDoublets) {
0122 for (unsigned int i = 0; i < layers.size(); i++) {
0123 for (unsigned int j = 0; j < 3; ++j) {
0124 auto vertexIndex = 0;
0125 auto foundVertex = std::find(g.theLayers.begin(), g.theLayers.end(), layers[i][j].name());
0126
0127 if (foundVertex == g.theLayers.end()) {
0128 vertexIndex = g.theLayers.size() - 1;
0129 } else {
0130 vertexIndex = foundVertex - g.theLayers.begin();
0131 }
0132
0133 if (j > 0) {
0134 auto innerVertex = std::find(g.theLayers.begin(), g.theLayers.end(), layers[i][j - 1].name());
0135
0136 CALayerPair tmpInnerLayerPair(innerVertex - g.theLayers.begin(), vertexIndex);
0137
0138 if (std::find(g.theLayerPairs.begin(), g.theLayerPairs.end(), tmpInnerLayerPair) == g.theLayerPairs.end()) {
0139 auto found = std::find_if(regionLayerPairs.begin(),
0140 regionLayerPairs.end(),
0141 [&](const IntermediateHitDoublets::LayerPairHitDoublets& pair) {
0142 return pair.innerLayerIndex() == layers[i][j - 1].index() &&
0143 pair.outerLayerIndex() == layers[i][j].index();
0144 });
0145 if (found != regionLayerPairs.end()) {
0146 hitDoublets.emplace_back(&(found->doublets()));
0147 g.theLayerPairs.push_back(tmpInnerLayerPair);
0148 g.theLayers[vertexIndex].theInnerLayers.push_back(innerVertex - g.theLayers.begin());
0149 innerVertex->theOuterLayers.push_back(vertexIndex);
0150 g.theLayers[vertexIndex].theInnerLayerPairs.push_back(g.theLayerPairs.size() - 1);
0151 innerVertex->theOuterLayerPairs.push_back(g.theLayerPairs.size() - 1);
0152 }
0153 }
0154 }
0155 }
0156 }
0157 }
0158 }
0159
0160 void CAHitTripletGenerator::hitNtuplets(const IntermediateHitDoublets& regionDoublets,
0161 std::vector<OrderedHitSeeds>& result,
0162 const SeedingLayerSetsHits& layers) {
0163 CAGraph g;
0164
0165 std::vector<const HitDoublets*> hitDoublets;
0166
0167 std::vector<CACell::CAntuplet> foundTriplets;
0168
0169 int index = 0;
0170 for (const auto& regionLayerPairs : regionDoublets) {
0171 const TrackingRegion& region = regionLayerPairs.region();
0172 hitDoublets.clear();
0173 foundTriplets.clear();
0174
0175 if (index == 0) {
0176 createGraphStructure(layers, g);
0177 caThetaCut.setCutValuesByLayerIds(g);
0178 caPhiCut.setCutValuesByLayerIds(g);
0179 } else {
0180 clearGraphStructure(layers, g);
0181 }
0182 fillGraph(layers, regionLayerPairs, g, hitDoublets);
0183 CellularAutomaton ca(g);
0184 ca.findTriplets(hitDoublets, foundTriplets, region, caThetaCut, caPhiCut, caHardPtCut);
0185
0186 auto& allCells = ca.getAllCells();
0187
0188 const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(*theField);
0189
0190
0191 std::array<float, 3> bc_r;
0192 std::array<float, 3> bc_z;
0193 std::array<float, 3> bc_errZ2;
0194 std::array<GlobalPoint, 3> gps;
0195 std::array<GlobalError, 3> ges;
0196 std::array<bool, 3> barrels;
0197
0198 unsigned int numberOfFoundTriplets = foundTriplets.size();
0199 for (unsigned int tripletId = 0; tripletId < numberOfFoundTriplets; ++tripletId) {
0200 OrderedHitTriplet tmpTriplet(allCells[foundTriplets[tripletId][0]].getInnerHit(),
0201 allCells[foundTriplets[tripletId][0]].getOuterHit(),
0202 allCells[foundTriplets[tripletId][1]].getOuterHit());
0203
0204 auto isBarrel = [](const unsigned id) -> bool { return id == PixelSubdetector::PixelBarrel; };
0205 for (unsigned int i = 0; i < 2; ++i) {
0206 auto const& ahit = allCells[foundTriplets[tripletId][i]].getInnerHit();
0207 gps[i] = ahit->globalPosition();
0208 ges[i] = ahit->globalPositionError();
0209 barrels[i] = isBarrel(ahit->geographicalId().subdetId());
0210 }
0211
0212 auto const& ahit = allCells[foundTriplets[tripletId][1]].getOuterHit();
0213 gps[2] = ahit->globalPosition();
0214 ges[2] = ahit->globalPositionError();
0215 barrels[2] = isBarrel(ahit->geographicalId().subdetId());
0216
0217 PixelRecoLineRZ line(gps[0], gps[2]);
0218 ThirdHitPredictionFromCircle predictionRPhi(gps[0], gps[2], extraHitRPhitolerance);
0219 const float curvature = predictionRPhi.curvature(ThirdHitPredictionFromCircle::Vector2D(gps[1].x(), gps[1].y()));
0220 const float abscurv = std::abs(curvature);
0221 const float thisMaxChi2 = maxChi2Eval.value(abscurv);
0222 float chi2 = std::numeric_limits<float>::quiet_NaN();
0223
0224
0225 if (useBendingCorrection) {
0226
0227 const float simpleCot = (gps.back().z() - gps.front().z()) / (gps.back().perp() - gps.front().perp());
0228 const float pt = 1.f / PixelRecoUtilities::inversePt(abscurv, *theField);
0229 for (int i = 0; i < 3; ++i) {
0230 const GlobalPoint& point = gps[i];
0231 const GlobalError& error = ges[i];
0232 bc_r[i] = sqrt(sqr(point.x() - region.origin().x()) + sqr(point.y() - region.origin().y()));
0233 bc_r[i] += pixelrecoutilities::LongitudinalBendingCorrection(pt, *theField)(bc_r[i]);
0234 bc_z[i] = point.z() - region.origin().z();
0235 bc_errZ2[i] = (barrels[i]) ? error.czz() : error.rerr(point) * sqr(simpleCot);
0236 }
0237 RZLine rzLine(bc_r, bc_z, bc_errZ2, RZLine::ErrZ2_tag());
0238 chi2 = rzLine.chi2();
0239 } else {
0240 RZLine rzLine(gps, ges, barrels);
0241 chi2 = rzLine.chi2();
0242 }
0243
0244 if (edm::isNotFinite(chi2) || chi2 > thisMaxChi2) {
0245 continue;
0246 }
0247
0248 if (theComparitor) {
0249 if (!theComparitor->compatible(tmpTriplet)) {
0250 continue;
0251 }
0252 }
0253 result[index].emplace_back(tmpTriplet);
0254 }
0255 index++;
0256 }
0257 }