PixelSeeding/src/CAHitQuadrupletGenerator.cc

0001 #include <functional>
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/CAHitQuadrupletGenerator.h"
0010 #include "RecoTracker/PixelSeeding/interface/ThirdHitPredictionFromCircle.h"
0011 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0012
0013 #include "RecoTracker/PixelSeeding/interface/CAGraph.h"
0014 #include "CellularAutomaton.h"
0015
0016 namespace {
0017   template <typename T>
0018   T sqr(T x) {
0019     return x * x;
0020   }
0021 }  // namespace
0022
0023 using namespace std;
0024
0025 constexpr unsigned int CAHitQuadrupletGenerator::minLayers;
0026
0027 CAHitQuadrupletGenerator::CAHitQuadrupletGenerator(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0028     : theFieldToken(iC.esConsumes()),
0029       extraHitRPhitolerance(cfg.getParameter<double>(
0030           "extraHitRPhitolerance")),  //extra window in ThirdHitPredictionFromCircle range (divide by R to get phi)
0031       maxChi2(cfg.getParameter<edm::ParameterSet>("maxChi2")),
0032       fitFastCircle(cfg.getParameter<bool>("fitFastCircle")),
0033       fitFastCircleChi2Cut(cfg.getParameter<bool>("fitFastCircleChi2Cut")),
0034       useBendingCorrection(cfg.getParameter<bool>("useBendingCorrection")),
0035       caThetaCut(cfg.getParameter<double>("CAThetaCut"),
0036                  cfg.getParameter<std::vector<edm::ParameterSet>>("CAThetaCut_byTriplets")),
0037       caPhiCut(cfg.getParameter<double>("CAPhiCut"),
0038                cfg.getParameter<std::vector<edm::ParameterSet>>("CAPhiCut_byTriplets")),
0039       caHardPtCut(cfg.getParameter<double>("CAHardPtCut")) {
0040   edm::ParameterSet comparitorPSet = cfg.getParameter<edm::ParameterSet>("SeedComparitorPSet");
0041   std::string comparitorName = comparitorPSet.getParameter<std::string>("ComponentName");
0042   if (comparitorName != "none") {
0043     theComparitor = SeedComparitorFactory::get()->create(comparitorName, comparitorPSet, iC);
0044   }
0045 }
0046
0047 void CAHitQuadrupletGenerator::fillDescriptions(edm::ParameterSetDescription& desc) {
0048   desc.add<double>("extraHitRPhitolerance", 0.1);
0049   desc.add<bool>("fitFastCircle", false);
0050   desc.add<bool>("fitFastCircleChi2Cut", false);
0051   desc.add<bool>("useBendingCorrection", false);
0052   desc.add<double>("CAThetaCut", 0.00125);
0053   desc.add<double>("CAPhiCut", 10);
0054
0055   edm::ParameterSetDescription validatorCACut;
0056   validatorCACut.add<string>("seedingLayers", "BPix1+BPix2+BPix3");
0057   validatorCACut.add<double>("cut", 0.00125);
0058   std::vector<edm::ParameterSet> defaultCACutVector;
0059   edm::ParameterSet defaultCACut;
0060   defaultCACut.addParameter<string>("seedingLayers", "");
0061   defaultCACut.addParameter<double>("cut", -1.);
0062   defaultCACutVector.push_back(defaultCACut);
0063   desc.addVPSet("CAThetaCut_byTriplets", validatorCACut, defaultCACutVector);
0064   desc.addVPSet("CAPhiCut_byTriplets", validatorCACut, defaultCACutVector);
0065
0066   desc.add<double>("CAHardPtCut", 0);
0067   desc.addOptional<bool>("CAOnlyOneLastHitPerLayerFilter")
0068       ->setComment(
0069           "Deprecated and has no effect. To be fully removed later when the parameter is no longer used in HLT "
0070           "configurations.");
0071   edm::ParameterSetDescription descMaxChi2;
0072   descMaxChi2.add<double>("pt1", 0.2);
0073   descMaxChi2.add<double>("pt2", 1.5);
0074   descMaxChi2.add<double>("value1", 500);
0075   descMaxChi2.add<double>("value2", 50);
0076   descMaxChi2.add<bool>("enabled", true);
0077   desc.add<edm::ParameterSetDescription>("maxChi2", descMaxChi2);
0078
0079   edm::ParameterSetDescription descComparitor;
0080   descComparitor.add<std::string>("ComponentName", "none");
0081   descComparitor.setAllowAnything();  // until we have moved SeedComparitor too to EDProducers
0082   desc.add<edm::ParameterSetDescription>("SeedComparitorPSet", descComparitor);
0083 }
0084
0085 void CAHitQuadrupletGenerator::initEvent(const edm::Event& ev, const edm::EventSetup& es) {
0086   if (theComparitor)
0087     theComparitor->init(ev, es);
0088   theField = &es.getData(theFieldToken);
0089 }
0090 namespace {
0091   void createGraphStructure(const SeedingLayerSetsHits& layers, CAGraph& g) {
0092     for (unsigned int i = 0; i < layers.size(); i++) {
0093       for (unsigned int j = 0; j < 4; ++j) {
0094         auto vertexIndex = 0;
0095         auto foundVertex = std::find(g.theLayers.begin(), g.theLayers.end(), layers[i][j].name());
0096         if (foundVertex == g.theLayers.end()) {
0097           g.theLayers.emplace_back(layers[i][j].name(), layers[i][j].detLayer()->seqNum(), layers[i][j].hits().size());
0098           vertexIndex = g.theLayers.size() - 1;
0099         } else {
0100           vertexIndex = foundVertex - g.theLayers.begin();
0101         }
0102         if (j == 0) {
0103           if (std::find(g.theRootLayers.begin(), g.theRootLayers.end(), vertexIndex) == g.theRootLayers.end()) {
0104             g.theRootLayers.emplace_back(vertexIndex);
0105           }
0106         }
0107       }
0108     }
0109   }
0110   void clearGraphStructure(const SeedingLayerSetsHits& layers, CAGraph& g) {
0111     g.theLayerPairs.clear();
0112     for (unsigned int i = 0; i < g.theLayers.size(); i++) {
0113       g.theLayers[i].theInnerLayers.clear();
0114       g.theLayers[i].theInnerLayerPairs.clear();
0115       g.theLayers[i].theOuterLayers.clear();
0116       g.theLayers[i].theOuterLayerPairs.clear();
0117       for (auto& v : g.theLayers[i].isOuterHitOfCell)
0118         v.clear();
0119     }
0120   }
0121   void fillGraph(const SeedingLayerSetsHits& layers,
0122                  const IntermediateHitDoublets::RegionLayerSets& regionLayerPairs,
0123                  CAGraph& g,
0124                  std::vector<const HitDoublets*>& hitDoublets) {
0125     for (unsigned int i = 0; i < layers.size(); i++) {
0126       for (unsigned int j = 0; j < 4; ++j) {
0127         auto vertexIndex = 0;
0128         auto foundVertex = std::find(g.theLayers.begin(), g.theLayers.end(), layers[i][j].name());
0129         if (foundVertex == g.theLayers.end()) {
0130           vertexIndex = g.theLayers.size() - 1;
0131         } else {
0132           vertexIndex = foundVertex - g.theLayers.begin();
0133         }
0134
0135         if (j > 0) {
0136           auto innerVertex = std::find(g.theLayers.begin(), g.theLayers.end(), layers[i][j - 1].name());
0137
0138           CALayerPair tmpInnerLayerPair(innerVertex - g.theLayers.begin(), vertexIndex);
0139
0140           if (std::find(g.theLayerPairs.begin(), g.theLayerPairs.end(), tmpInnerLayerPair) == g.theLayerPairs.end()) {
0141             auto found = std::find_if(regionLayerPairs.begin(),
0142                                       regionLayerPairs.end(),
0143                                       [&](const IntermediateHitDoublets::LayerPairHitDoublets& pair) {
0144                                         return pair.innerLayerIndex() == layers[i][j - 1].index() &&
0145                                                pair.outerLayerIndex() == layers[i][j].index();
0146                                       });
0147             if (found != regionLayerPairs.end()) {
0148               hitDoublets.emplace_back(&(found->doublets()));
0149               g.theLayerPairs.push_back(tmpInnerLayerPair);
0150               g.theLayers[vertexIndex].theInnerLayers.push_back(innerVertex - g.theLayers.begin());
0151               innerVertex->theOuterLayers.push_back(vertexIndex);
0152               g.theLayers[vertexIndex].theInnerLayerPairs.push_back(g.theLayerPairs.size() - 1);
0153               innerVertex->theOuterLayerPairs.push_back(g.theLayerPairs.size() - 1);
0154             }
0155           }
0156         }
0157       }
0158     }
0159   }
0160 }  // namespace
0161
0162 void CAHitQuadrupletGenerator::hitNtuplets(const IntermediateHitDoublets& regionDoublets,
0163                                            std::vector<OrderedHitSeeds>& result,
0164                                            const SeedingLayerSetsHits& layers) {
0165   CAGraph g;
0166
0167   std::vector<const HitDoublets*> hitDoublets;
0168
0169   const int numberOfHitsInNtuplet = 4;
0170   std::vector<CACell::CAntuplet> foundQuadruplets;
0171
0172   int index = 0;
0173   for (const auto& regionLayerPairs : regionDoublets) {
0174     const TrackingRegion& region = regionLayerPairs.region();
0175     hitDoublets.clear();
0176     foundQuadruplets.clear();
0177     if (index == 0) {
0178       createGraphStructure(layers, g);
0179       caThetaCut.setCutValuesByLayerIds(g);
0180       caPhiCut.setCutValuesByLayerIds(g);
0181     } else {
0182       clearGraphStructure(layers, g);
0183     }
0184
0185     fillGraph(layers, regionLayerPairs, g, hitDoublets);
0186
0187     CellularAutomaton ca(g);
0188
0189     ca.createAndConnectCells(hitDoublets, region, caThetaCut, caPhiCut, caHardPtCut);
0190
0191     ca.evolve(numberOfHitsInNtuplet);
0192
0193     ca.findNtuplets(foundQuadruplets, numberOfHitsInNtuplet);
0194
0195     auto& allCells = ca.getAllCells();
0196
0197     const QuantityDependsPtEval maxChi2Eval = maxChi2.evaluator(*theField);
0198
0199     // re-used thoughout
0200     std::array<float, 4> bc_r;
0201     std::array<float, 4> bc_z;
0202     std::array<float, 4> bc_errZ2;
0203     std::array<GlobalPoint, 4> gps;
0204     std::array<GlobalError, 4> ges;
0205     std::array<bool, 4> barrels;
0206
0207     unsigned int numberOfFoundQuadruplets = foundQuadruplets.size();
0208
0209     // Loop over quadruplets
0210     for (unsigned int quadId = 0; quadId < numberOfFoundQuadruplets; ++quadId) {
0211       auto isBarrel = [](const unsigned id) -> bool { return id == PixelSubdetector::PixelBarrel; };
0212       for (unsigned int i = 0; i < 3; ++i) {
0213         auto const& ahit = allCells[foundQuadruplets[quadId][i]].getInnerHit();
0214         gps[i] = ahit->globalPosition();
0215         ges[i] = ahit->globalPositionError();
0216         barrels[i] = isBarrel(ahit->geographicalId().subdetId());
0217       }
0218
0219       auto const& ahit = allCells[foundQuadruplets[quadId][2]].getOuterHit();
0220       gps[3] = ahit->globalPosition();
0221       ges[3] = ahit->globalPositionError();
0222       barrels[3] = isBarrel(ahit->geographicalId().subdetId());
0223       // TODO:
0224       // - if we decide to always do the circle fit for 4 hits, we don't
0225       //   need ThirdHitPredictionFromCircle for the curvature; then we
0226       //   could remove extraHitRPhitolerance configuration parameter
0227       ThirdHitPredictionFromCircle predictionRPhi(gps[0], gps[2], extraHitRPhitolerance);
0228       const float curvature = predictionRPhi.curvature(ThirdHitPredictionFromCircle::Vector2D(gps[1].x(), gps[1].y()));
0229       const float abscurv = std::abs(curvature);
0230       const float thisMaxChi2 = maxChi2Eval.value(abscurv);
0231       if (theComparitor) {
0232         SeedingHitSet tmpTriplet(allCells[foundQuadruplets[quadId][0]].getInnerHit(),
0233                                  allCells[foundQuadruplets[quadId][2]].getInnerHit(),
0234                                  allCells[foundQuadruplets[quadId][2]].getOuterHit());
0235
0236         if (!theComparitor->compatible(tmpTriplet)) {
0237           continue;
0238         }
0239       }
0240
0241       float chi2 = std::numeric_limits<float>::quiet_NaN();
0242       // TODO: Do we have any use case to not use bending correction?
0243       if (useBendingCorrection) {
0244         // Following PixelFitterByConformalMappingAndLine
0245         const float simpleCot = (gps.back().z() - gps.front().z()) / (gps.back().perp() - gps.front().perp());
0246         const float pt = 1.f / PixelRecoUtilities::inversePt(abscurv, *theField);
0247         for (int i = 0; i < 4; ++i) {
0248           const GlobalPoint& point = gps[i];
0249           const GlobalError& error = ges[i];
0250           bc_r[i] = sqrt(sqr(point.x() - region.origin().x()) + sqr(point.y() - region.origin().y()));
0251           bc_r[i] += pixelrecoutilities::LongitudinalBendingCorrection(pt, *theField)(bc_r[i]);
0252           bc_z[i] = point.z() - region.origin().z();
0253           bc_errZ2[i] = (barrels[i]) ? error.czz() : error.rerr(point) * sqr(simpleCot);
0254         }
0255         RZLine rzLine(bc_r, bc_z, bc_errZ2, RZLine::ErrZ2_tag());
0256         chi2 = rzLine.chi2();
0257       } else {
0258         RZLine rzLine(gps, ges, barrels);
0259         chi2 = rzLine.chi2();
0260       }
0261       if (edm::isNotFinite(chi2) || chi2 > thisMaxChi2) {
0262         continue;
0263       }
0264       // TODO: Do we have any use case to not use circle fit? Maybe
0265       // HLT where low-pT inefficiency is not a problem?
0266       if (fitFastCircle) {
0267         FastCircleFit c(gps, ges);
0268         chi2 += c.chi2();
0269         if (edm::isNotFinite(chi2))
0270           continue;
0271         if (fitFastCircleChi2Cut && chi2 > thisMaxChi2)
0272           continue;
0273       }
0274       result[index].emplace_back(allCells[foundQuadruplets[quadId][0]].getInnerHit(),
0275                                  allCells[foundQuadruplets[quadId][1]].getInnerHit(),
0276                                  allCells[foundQuadruplets[quadId][2]].getInnerHit(),
0277                                  allCells[foundQuadruplets[quadId][2]].getOuterHit());
0278     }
0279     index++;
0280   }
0281 }