File indexing completed on 2024-04-06 12:28:34
0001 #include "PixelTripletNoTipGenerator.h"
0002 #include "FWCore/Framework/interface/ConsumesCollector.h"
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005
0006 #include "RecoTracker/TkHitPairs/interface/RecHitsSortedInPhi.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "RecoTracker/TkMSParametrization/interface/PixelRecoUtilities.h"
0009 #include "ThirdHitPredictionFromInvLine.h"
0010 #include "ThirdHitZPrediction.h"
0011
0012 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisation.h"
0013 #include "RecoTracker/TkHitPairs/interface/HitPairGeneratorFromLayerPair.h"
0014
0015 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0016 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0017 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0018 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
0019
0020 typedef PixelRecoRange<float> Range;
0021 template <class T>
0022 T sqr(T t) {
0023 return t * t;
0024 }
0025
0026 using namespace std;
0027
0028 PixelTripletNoTipGenerator::PixelTripletNoTipGenerator(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0029 : HitTripletGeneratorFromPairAndLayers(cfg),
0030 fieldToken_(iC.esConsumes()),
0031 msmakerToken_(iC.esConsumes()),
0032 extraHitRZtolerance(cfg.getParameter<double>("extraHitRZtolerance")),
0033 extraHitRPhitolerance(cfg.getParameter<double>("extraHitRPhitolerance")),
0034 extraHitPhiToleranceForPreFiltering(cfg.getParameter<double>("extraHitPhiToleranceForPreFiltering")),
0035 theNSigma(cfg.getParameter<double>("nSigma")),
0036 theChi2Cut(cfg.getParameter<double>("chi2Cut")) {}
0037
0038 PixelTripletNoTipGenerator::~PixelTripletNoTipGenerator() {}
0039
0040 void PixelTripletNoTipGenerator::hitTriplets(const TrackingRegion& region,
0041 OrderedHitTriplets& result,
0042 const edm::Event& ev,
0043 const edm::EventSetup& es,
0044 const SeedingLayerSetsHits::SeedingLayerSet& pairLayers,
0045 const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers) {
0046
0047
0048
0049
0050
0051
0052
0053
0054 const GlobalPoint& bsPoint = region.origin();
0055 double errorXY = region.originRBound();
0056 GlobalVector shift = bsPoint - GlobalPoint(0., 0., 0.);
0057
0058 OrderedHitPairs pairs;
0059 pairs.reserve(30000);
0060 OrderedHitPairs::const_iterator ip;
0061 thePairGenerator->hitPairs(region, pairs, ev, es, pairLayers);
0062
0063 if (pairs.empty())
0064 return;
0065
0066 const DetLayer* firstLayer = thePairGenerator->innerLayer(pairLayers).detLayer();
0067 const DetLayer* secondLayer = thePairGenerator->outerLayer(pairLayers).detLayer();
0068 if (!firstLayer || !secondLayer)
0069 return;
0070
0071 int size = thirdLayers.size();
0072
0073 const RecHitsSortedInPhi** thirdHitMap = new const RecHitsSortedInPhi*[size];
0074 for (int il = 0; il <= size - 1; il++) {
0075 thirdHitMap[il] = &(*theLayerCache)(thirdLayers[il], region);
0076 }
0077
0078 const auto& field = es.getData(fieldToken_);
0079 const auto& msmaker = es.getData(msmakerToken_);
0080
0081 MultipleScatteringParametrisation sigma1RPhi = msmaker.parametrisation(firstLayer);
0082 MultipleScatteringParametrisation sigma2RPhi = msmaker.parametrisation(secondLayer);
0083
0084 typedef RecHitsSortedInPhi::Hit Hit;
0085 for (ip = pairs.begin(); ip != pairs.end(); ip++) {
0086 GlobalPoint p1((*ip).inner()->globalPosition() - shift);
0087 GlobalPoint p2((*ip).outer()->globalPosition() - shift);
0088
0089 ThirdHitPredictionFromInvLine predictionRPhiTMP(p1, p2);
0090 double pt_p1p2 = 1. / PixelRecoUtilities::inversePt(predictionRPhiTMP.curvature(), field);
0091
0092 PixelRecoPointRZ point1(p1.perp(), p1.z());
0093 PixelRecoPointRZ point2(p2.perp(), p2.z());
0094
0095 PixelRecoLineRZ line(point1, point2);
0096 double msRPhi1 = sigma1RPhi(pt_p1p2, line.cotLine(), 0.f);
0097 double msRPhi2 = sigma2RPhi(pt_p1p2, line.cotLine(), point1);
0098 double sinTheta = 1 / sqrt(1 + sqr(line.cotLine()));
0099 double cosTheta = fabs(line.cotLine()) / sqrt(1 + sqr(line.cotLine()));
0100
0101 double p1_errorRPhi = sqrt(sqr((*ip).inner()->errorGlobalRPhi()) + sqr(msRPhi1) + sqr(errorXY));
0102 double p2_errorRPhi = sqrt(sqr((*ip).outer()->errorGlobalRPhi()) + sqr(msRPhi2) + sqr(errorXY));
0103
0104 ThirdHitPredictionFromInvLine predictionRPhi(p1, p2, p1_errorRPhi, p2_errorRPhi);
0105
0106 for (int il = 0; il <= size - 1; il++) {
0107 const DetLayer* layer = thirdLayers[il].detLayer();
0108 bool barrelLayer = (layer->location() == GeomDetEnumerators::barrel);
0109 MultipleScatteringParametrisation sigma3RPhi = msmaker.parametrisation(layer);
0110 double msRPhi3 = sigma3RPhi(pt_p1p2, line.cotLine(), point2);
0111
0112 Range rRange;
0113 if (barrelLayer) {
0114 const BarrelDetLayer& bl = dynamic_cast<const BarrelDetLayer&>(*layer);
0115 float halfThickness = bl.surface().bounds().thickness() / 2;
0116 float radius = bl.specificSurface().radius();
0117 rRange = Range(radius - halfThickness, radius + halfThickness);
0118 } else {
0119 const ForwardDetLayer& fl = dynamic_cast<const ForwardDetLayer&>(*layer);
0120 float halfThickness = fl.surface().bounds().thickness() / 2;
0121 float zLayer = fl.position().z();
0122 float zMin = zLayer - halfThickness;
0123 float zMax = zLayer + halfThickness;
0124 GlobalVector dr = p2 - p1;
0125 GlobalPoint p3_a = p2 + dr * (zMin - p2.z()) / dr.z();
0126 GlobalPoint p3_b = p2 + dr * (zMax - p2.z()) / dr.z();
0127 if (zLayer * p3_a.z() < 0)
0128 continue;
0129 rRange = Range(p3_a.perp(), p3_b.perp());
0130 rRange.sort();
0131 }
0132 double displacment = shift.perp();
0133 GlobalPoint crossing1 = predictionRPhi.crossing(rRange.min() - displacment) + shift;
0134 GlobalPoint crossing2 = predictionRPhi.crossing(rRange.max() + displacment) + shift;
0135 float c1_phi = crossing1.phi();
0136 float c2_phi = crossing2.phi();
0137 if (c2_phi < c1_phi)
0138 swap(c1_phi, c2_phi);
0139 if (c2_phi - c1_phi > M_PI) {
0140 c2_phi -= 2 * M_PI;
0141 swap(c1_phi, c2_phi);
0142 }
0143 double extraAngle = (displacment + theNSigma * msRPhi3) / rRange.min() + extraHitPhiToleranceForPreFiltering;
0144 c1_phi -= extraAngle;
0145 c2_phi += extraAngle;
0146 vector<Hit> thirdHits = thirdHitMap[il]->hits(c1_phi, c2_phi);
0147
0148 typedef vector<Hit>::const_iterator IH;
0149 for (IH th = thirdHits.begin(), eh = thirdHits.end(); th < eh; ++th) {
0150 GlobalPoint p3((*th)->globalPosition() - shift);
0151 double p3_errorRPhi = sqrt(sqr((*th)->errorGlobalRPhi()) + sqr(msRPhi3) + sqr(errorXY));
0152
0153 predictionRPhi.add(p3, p3_errorRPhi);
0154
0155 double curvature = predictionRPhi.curvature();
0156
0157 ThirdHitZPrediction zPrediction((*ip).inner()->globalPosition(),
0158 sqrt(sqr((*ip).inner()->errorGlobalR()) + sqr(msRPhi1 / cosTheta)),
0159 sqrt(sqr((*ip).inner()->errorGlobalZ()) + sqr(msRPhi1 / sinTheta)),
0160 (*ip).outer()->globalPosition(),
0161 sqrt(sqr((*ip).outer()->errorGlobalR()) + sqr(msRPhi2 / cosTheta)),
0162 sqrt(sqr((*ip).outer()->errorGlobalZ()) + sqr(msRPhi2 / sinTheta)),
0163 curvature,
0164 theNSigma);
0165 ThirdHitZPrediction::Range zRange =
0166 zPrediction((*th)->globalPosition(), sqrt(sqr((*th)->errorGlobalR())) + sqr(msRPhi3 / cosTheta));
0167
0168 double z3Hit = (*th)->globalPosition().z();
0169 double z3HitError =
0170 theNSigma * (sqrt(sqr((*th)->errorGlobalZ()) + sqr(msRPhi3 / sinTheta))) + extraHitRZtolerance;
0171 Range hitZRange(z3Hit - z3HitError, z3Hit + z3HitError);
0172 bool inside = hitZRange.hasIntersection(zRange);
0173
0174 double curvatureMS = PixelRecoUtilities::curvature(1. / region.ptMin(), field);
0175 bool ptCut = (predictionRPhi.curvature() - theNSigma * predictionRPhi.errorCurvature() < curvatureMS);
0176 bool chi2Cut = (predictionRPhi.chi2() < theChi2Cut);
0177 if (inside && ptCut && chi2Cut) {
0178 if (theMaxElement != 0 && result.size() >= theMaxElement) {
0179 result.clear();
0180 edm::LogError("TooManyTriplets") << " number of triples exceed maximum. no triplets produced.";
0181 delete[] thirdHitMap;
0182 return;
0183 }
0184 result.push_back(OrderedHitTriplet((*ip).inner(), (*ip).outer(), *th));
0185 }
0186 predictionRPhi.remove(p3, p3_errorRPhi);
0187 }
0188 }
0189 }
0190 delete[] thirdHitMap;
0191 }
0192 void PixelTripletNoTipGenerator::hitTriplets(const TrackingRegion& region,
0193 OrderedHitTriplets& result,
0194 const edm::EventSetup& es,
0195 const HitDoublets& doublets,
0196 const RecHitsSortedInPhi** thirdHitMap,
0197 const std::vector<const DetLayer*>& thirdLayerDetLayer,
0198 const int nThirdLayers) {
0199 throw cms::Exception("Error") << "PixelTripletNoTipGenerator::hitTriplets is not implemented \n";
0200 }