PixelSeeding/src/CAHitTripletGenerator.cc

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 }  // namespace
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")),  //extra window in ThirdHitPredictionFromCircle range (divide by R to get phi)
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();  // until we have moved SeedComparitor too to EDProducers
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 }  // namespace
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     // re-used thoughout, need to be vectors because of RZLine interface
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       // TODO: Do we have any use case to not use bending correction?
0224
0225       if (useBendingCorrection) {
0226         // Following PixelFitterByConformalMappingAndLine
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 }