PixelLowPtUtilities/src/PixelTripletLowPtGenerator.cc

0001 #include "RecoPixelVertexing/PixelLowPtUtilities/interface/PixelTripletLowPtGenerator.h"
0002 #include "RecoPixelVertexing/PixelLowPtUtilities/interface/ThirdHitPrediction.h"
0003 #include "RecoPixelVertexing/PixelLowPtUtilities/interface/TripletFilter.h"
0004 #include "RecoPixelVertexing/PixelLowPtUtilities/interface/HitInfo.h"
0005
0006 #include "RecoTracker/TkMSParametrization/interface/PixelRecoPointRZ.h"
0007 #include "RecoTracker/TkHitPairs/interface/HitPairGeneratorFromLayerPair.h"
0008
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "DataFormats/SiPixelCluster/interface/SiPixelClusterShapeCache.h"
0013
0014 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0015 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.h"
0016 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0017 #include "RecoTracker/Record/interface/TrackerMultipleScatteringRecord.h"
0018
0019 #undef Debug
0020
0021 using namespace std;
0022
0023 /*****************************************************************************/
0024 PixelTripletLowPtGenerator::PixelTripletLowPtGenerator(const edm::ParameterSet& cfg, edm::ConsumesCollector& iC)
0025     : HitTripletGeneratorFromPairAndLayers(),  // no theMaxElement used in this class
0026       m_geomToken(iC.esConsumes<TrackerGeometry, TrackerDigiGeometryRecord>()),
0027       m_topoToken(iC.esConsumes<TrackerTopology, TrackerTopologyRcd>()),
0028       m_magfieldToken(iC.esConsumes()),
0029       m_ttrhBuilderToken(iC.esConsumes(edm::ESInputTag("", cfg.getParameter<string>("TTRHBuilder")))),
0030       m_msmakerToken(iC.esConsumes()),
0031       theTracker(nullptr),
0032       theClusterShapeCacheToken(
0033           iC.consumes<SiPixelClusterShapeCache>(cfg.getParameter<edm::InputTag>("clusterShapeCacheSrc"))) {
0034   checkMultipleScattering = cfg.getParameter<bool>("checkMultipleScattering");
0035   nSigMultipleScattering = cfg.getParameter<double>("nSigMultipleScattering");
0036   checkClusterShape = cfg.getParameter<bool>("checkClusterShape");
0037   rzTolerance = cfg.getParameter<double>("rzTolerance");
0038   maxAngleRatio = cfg.getParameter<double>("maxAngleRatio");
0039 }
0040
0041 /*****************************************************************************/
0042 PixelTripletLowPtGenerator::~PixelTripletLowPtGenerator() {}
0043
0044 /*****************************************************************************/
0045
0046 /*****************************************************************************/
0047 void PixelTripletLowPtGenerator::getTracker(const edm::EventSetup& es) {
0048   if (theTracker == nullptr) {
0049     // Get tracker geometry
0050     theTracker = &es.getData(m_geomToken);
0051   }
0052
0053   if (!theFilter) {
0054     theFilter = std::make_unique<TripletFilter>(es);
0055   }
0056 }
0057
0058 /*****************************************************************************/
0059 GlobalPoint PixelTripletLowPtGenerator::getGlobalPosition(const TrackingRecHit* recHit) {
0060   DetId detId = recHit->geographicalId();
0061
0062   return theTracker->idToDet(detId)->toGlobal(recHit->localPosition());
0063 }
0064
0065 /*****************************************************************************/
0066 void PixelTripletLowPtGenerator::hitTriplets(const TrackingRegion& region,
0067                                              OrderedHitTriplets& result,
0068                                              const edm::Event& ev,
0069                                              const edm::EventSetup& es,
0070                                              const SeedingLayerSetsHits::SeedingLayerSet& pairLayers,
0071                                              const std::vector<SeedingLayerSetsHits::SeedingLayer>& thirdLayers) {
0072   //Retrieve tracker topology from geometry
0073   const TrackerTopology* tTopo = &es.getData(m_topoToken);
0074   const auto& magfield = es.getData(m_magfieldToken);
0075   const auto& ttrhBuilder = es.getData(m_ttrhBuilderToken);
0076   const auto& msmaker = es.getData(m_msmakerToken);
0077
0078   edm::Handle<SiPixelClusterShapeCache> clusterShapeCache;
0079   ev.getByToken(theClusterShapeCacheToken, clusterShapeCache);
0080
0081   // Generate pairs
0082   OrderedHitPairs pairs;
0083   pairs.reserve(30000);
0084   thePairGenerator->hitPairs(region, pairs, ev, es, pairLayers);
0085
0086   if (pairs.empty())
0087     return;
0088
0089   int size = thirdLayers.size();
0090
0091   // Set aliases
0092   const RecHitsSortedInPhi** thirdHitMap = new const RecHitsSortedInPhi*[size];
0093   for (int il = 0; il < size; il++)
0094     thirdHitMap[il] = &(*theLayerCache)(thirdLayers[il], region);
0095
0096   // Get tracker
0097   getTracker(es);
0098
0099   // Look at all generated pairs
0100   for (OrderedHitPairs::const_iterator ip = pairs.begin(); ip != pairs.end(); ip++) {
0101     // Fill rechits and points
0102     vector<const TrackingRecHit*> recHits(3);
0103     vector<GlobalPoint> points(3);
0104
0105     recHits[0] = (*ip).inner()->hit();
0106     recHits[1] = (*ip).outer()->hit();
0107
0108 #ifdef Debug
0109     cerr << " RecHits " + HitInfo::getInfo(*recHits[0]) + HitInfo::getInfo(*recHits[1]) << endl;
0110 #endif
0111
0112     for (int i = 0; i < 2; i++)
0113       points[i] = getGlobalPosition(recHits[i]);
0114
0115     // Initialize helix prediction
0116     ThirdHitPrediction thePrediction(
0117         region, points[0], points[1], magfield, ttrhBuilder, nSigMultipleScattering, maxAngleRatio);
0118
0119     // Look at all layers
0120     for (int il = 0; il < size; il++) {
0121       const DetLayer* layer = thirdLayers[il].detLayer();
0122
0123 #ifdef Debug
0124       cerr << "  check layer " << layer->subDetector() << " " << layer->location() << endl;
0125 #endif
0126
0127       // Get ranges for the third hit
0128       float phi[2], rz[2];
0129       thePrediction.getRanges(layer, phi, rz);
0130
0131       PixelRecoRange<float> phiRange(phi[0], phi[1]);
0132       PixelRecoRange<float> rzRange(rz[0] - rzTolerance, rz[1] + rzTolerance);
0133
0134       // Get third hit candidates from cache
0135       typedef RecHitsSortedInPhi::Hit Hit;
0136       vector<Hit> thirdHits = thirdHitMap[il]->hits(phiRange.min(), phiRange.max());
0137       typedef vector<Hit>::const_iterator IH;
0138
0139       for (IH th = thirdHits.begin(), eh = thirdHits.end(); th < eh; ++th) {
0140         // Fill rechit and point
0141         recHits[2] = (*th)->hit();
0142         points[2] = getGlobalPosition(recHits[2]);
0143
0144 #ifdef Debug
0145         cerr << "  third hit " + HitInfo::getInfo(*recHits[2]) << endl;
0146 #endif
0147
0148         // Check if third hit is compatible with multiple scattering
0149         vector<GlobalVector> globalDirs;
0150         if (thePrediction.isCompatibleWithMultipleScattering(points[2], recHits, globalDirs, msmaker) == false) {
0151 #ifdef Debug
0152           cerr << "  not compatible: multiple scattering" << endl;
0153 #endif
0154           if (checkMultipleScattering)
0155             continue;
0156         }
0157
0158         // Convert to localDirs
0159         /*
0160         vector<LocalVector> localDirs;
0161         vector<GlobalVector>::const_iterator globalDir = globalDirs.begin();
0162         for(vector<const TrackingRecHit *>::const_iterator
0163                                             recHit  = recHits.begin();
0164                                             recHit != recHits.end(); recHit++)
0165         {
0166           localDirs.push_back(theTracker->idToDet(
0167                              (*recHit)->geographicalId())->toLocal(*globalDir));
0168           globalDir++;
0169         }
0170 */
0171
0172         // Check if the cluster shapes are compatible with thrusts
0173         if (checkClusterShape) {
0174           if (!theFilter->checkTrack(recHits, globalDirs, tTopo, *clusterShapeCache)) {
0175 #ifdef Debug
0176             cerr << "  not compatible: cluster shape" << endl;
0177 #endif
0178             continue;
0179           }
0180         }
0181
0182         // All checks passed, put triplet back
0183         result.push_back(OrderedHitTriplet((*ip).inner(), (*ip).outer(), *th));
0184       }
0185     }
0186   }
0187   delete[] thirdHitMap;
0188
0189   return;
0190 }
0191 void PixelTripletLowPtGenerator::hitTriplets(const TrackingRegion& region,
0192                                              OrderedHitTriplets& result,
0193                                              const edm::EventSetup& es,
0194                                              const HitDoublets& doublets,
0195                                              const RecHitsSortedInPhi** thirdHitMap,
0196                                              const std::vector<const DetLayer*>& thirdLayerDetLayer,
0197                                              const int nThirdLayers) {
0198   throw cms::Exception("Error") << "PixelTripletLowPtGenerator::hitTriplets is not implemented \n";
0199 }