HiTracking/plugins/HITrackingRegionProducer.h

0001 #ifndef RecoHI_HiTracking_HITrackingRegionProducer_H
0002 #define RecoHI_HiTracking_HITrackingRegionProducer_H
0003
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/Framework/interface/ConsumesCollector.h"
0006
0007 #include "RecoTracker/TkTrackingRegions/interface/TrackingRegionProducer.h"
0008 #include "RecoTracker/TkTrackingRegions/interface/GlobalTrackingRegion.h"
0009 #include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
0010 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "DataFormats/Common/interface/DetSetVector.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0014 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0015 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0016 #include "MagneticField/Engine/interface/MagneticField.h"
0017 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0018 #include "RecoTracker/Record/interface/TrackerMultipleScatteringRecord.h"
0019 #include "RecoTracker/TkMSParametrization/interface/MultipleScatteringParametrisationMaker.h"
0020
0021 #include "TMath.h"
0022
0023 class HITrackingRegionProducer : public TrackingRegionProducer {
0024 public:
0025   HITrackingRegionProducer(const edm::ParameterSet& cfg, edm::ConsumesCollector&& iC)
0026       : theTtopoToken(iC.esConsumes()), theFieldToken(iC.esConsumes()) {
0027     edm::ParameterSet regionPSet = cfg.getParameter<edm::ParameterSet>("RegionPSet");
0028
0029     thePtMin = regionPSet.getParameter<double>("ptMin");
0030     theOriginRadius = regionPSet.getParameter<double>("originRadius");
0031     theOriginHalfLength = regionPSet.getParameter<double>("originHalfLength");
0032     double xPos = regionPSet.getParameter<double>("originXPos");
0033     double yPos = regionPSet.getParameter<double>("originYPos");
0034     double zPos = regionPSet.getParameter<double>("originZPos");
0035     double xDir = regionPSet.getParameter<double>("directionXCoord");
0036     double yDir = regionPSet.getParameter<double>("directionYCoord");
0037     double zDir = regionPSet.getParameter<double>("directionZCoord");
0038     thePrecise = regionPSet.getParameter<bool>("precise");
0039     theSiPixelRecHits = regionPSet.getParameter<edm::InputTag>("siPixelRecHits");
0040     theSiPixelRecHitsToken = iC.consumes<SiPixelRecHitCollection>(theSiPixelRecHits);
0041     theOrigin = GlobalPoint(xPos, yPos, zPos);
0042     theDirection = GlobalVector(xDir, yDir, zDir);
0043     if (thePrecise) {
0044       theMSMakerToken = iC.esConsumes();
0045     }
0046   }
0047
0048   ~HITrackingRegionProducer() override = default;
0049
0050   int estimateMultiplicity(const edm::Event& ev, const edm::EventSetup& es) const {
0051     //rechits
0052     edm::Handle<SiPixelRecHitCollection> recHitColl;
0053     ev.getByToken(theSiPixelRecHitsToken, recHitColl);
0054
0055     //Retrieve tracker topology from geometry
0056     const auto& tTopo = es.getData(theTtopoToken);
0057
0058     int numRecHits = 0;
0059     //FIXME: this can be optimized quite a bit by looping only on the per-det 'items' of DetSetVector
0060     for (SiPixelRecHitCollection::const_iterator recHitIdIterator = recHitColl->begin(),
0061                                                  recHitIdIteratorEnd = recHitColl->end();
0062          recHitIdIterator != recHitIdIteratorEnd;
0063          recHitIdIterator++) {
0064       SiPixelRecHitCollection::DetSet hits = *recHitIdIterator;
0065       DetId detId = DetId(hits.detId());      // Get the Detid object
0066       unsigned int detType = detId.det();     // det type, tracker=1
0067       unsigned int subid = detId.subdetId();  //subdetector type, barrel=1, fpix=2
0068
0069       unsigned int layer = 0;
0070       layer = tTopo.pxbLayer(detId);
0071       if (detType == 1 && subid == 1 && layer == 1) {
0072         numRecHits += hits.size();
0073       }
0074     }
0075     return numRecHits;
0076   }
0077
0078   std::vector<std::unique_ptr<TrackingRegion> > regions(const edm::Event& ev,
0079                                                         const edm::EventSetup& es) const override {
0080     int estMult = estimateMultiplicity(ev, es);
0081
0082     // fit from MC information
0083     float aa = 1.90935e-04;
0084     float bb = -2.90167e-01;
0085     float cc = 3.86125e+02;
0086
0087     float estTracks = aa * estMult * estMult + bb * estMult + cc;
0088
0089     LogTrace("heavyIonHLTVertexing") << "[HIVertexing]";
0090     LogTrace("heavyIonHLTVertexing") << " [HIVertexing: hits in the 1. layer:" << estMult << "]";
0091     LogTrace("heavyIonHLTVertexing") << " [HIVertexing: estimated number of tracks:" << estTracks << "]";
0092
0093     float regTracking = 400.;  //if we have more tracks -> regional tracking
0094     float etaB = 10.;
0095     float phiB = TMath::Pi() / 2.;
0096
0097     float decEta = estTracks / 600.;
0098     etaB = 2.5 / decEta;
0099
0100     if (estTracks > regTracking) {
0101       LogTrace("heavyIonHLTVertexing") << " [HIVertexing: Regional Tracking]";
0102       LogTrace("heavyIonHLTVertexing") << "  [Regional Tracking: eta range: -" << etaB << ", " << etaB << "]";
0103       LogTrace("heavyIonHLTVertexing") << "  [Regional Tracking: phi range: -" << phiB << ", " << phiB << "]";
0104       LogTrace("heavyIonHLTVertexing") << "  [Regional Tracking: factor of decrease: " << decEta * 2.
0105                                        << "]";  // 2:from phi
0106     }
0107
0108     // tracking region selection
0109     std::vector<std::unique_ptr<TrackingRegion> > result;
0110     if (estTracks > regTracking) {  // regional tracking
0111       const auto& field = es.getData(theFieldToken);
0112       const MultipleScatteringParametrisationMaker* msmaker = nullptr;
0113       if (thePrecise) {
0114         msmaker = &es.getData(theMSMakerToken);
0115       }
0116
0117       result.push_back(std::make_unique<RectangularEtaPhiTrackingRegion>(theDirection,
0118                                                                          theOrigin,
0119                                                                          thePtMin,
0120                                                                          theOriginRadius,
0121                                                                          theOriginHalfLength,
0122                                                                          etaB,
0123                                                                          phiB,
0124                                                                          field,
0125                                                                          msmaker,
0126                                                                          thePrecise));
0127     } else {  // global tracking
0128       LogTrace("heavyIonHLTVertexing") << " [HIVertexing: Global Tracking]";
0129       result.push_back(std::make_unique<GlobalTrackingRegion>(
0130           thePtMin, theOrigin, theOriginRadius, theOriginHalfLength, thePrecise));
0131     }
0132     return result;
0133   }
0134
0135 private:
0136   edm::InputTag theSiPixelRecHits;
0137   edm::EDGetTokenT<SiPixelRecHitCollection> theSiPixelRecHitsToken;
0138   edm::ESGetToken<TrackerTopology, IdealGeometryRecord> theTtopoToken;
0139   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> theFieldToken;
0140   edm::ESGetToken<MultipleScatteringParametrisationMaker, TrackerMultipleScatteringRecord> theMSMakerToken;
0141   double thePtMin;
0142   GlobalPoint theOrigin;
0143   double theOriginRadius;
0144   double theOriginHalfLength;
0145   bool thePrecise;
0146   GlobalVector theDirection;
0147 };
0148
0149 #endif