PixelTriplets/plugins/PixelTripletNoTipGenerator.cc

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   //edm::Handle<reco::BeamSpot> bsHandle;
0048   //ev.getByLabel( theBeamSpotTag, bsHandle);
0049   //if(!bsHandle.isValid()) return;
0050   //const reco::BeamSpot & bs = *bsHandle;
0051   //double errorXY = sqrt( sqr(bs.BeamWidthX()) + sqr(bs.BeamWidthY()) );
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 }