MeasurementDet/test/MeasurementTrackerTest.cc

0001 #include "DataFormats/Common/interface/Handle.h"
0002 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0003 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0004 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0005 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0006 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0007 #include "FWCore/Framework/interface/Event.h"
0008 #include "FWCore/Framework/interface/EventSetup.h"
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Framework/interface/MakerMacros.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0013 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0014 #include "MagneticField/Engine/interface/MagneticField.h"
0015 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0016 #include "MagneticField/VolumeGeometry/interface/MagVolumeOutsideValidity.h"
0017 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0018 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0019 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
0020 #include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h"
0021 #include "TrackingTools/DetLayers/interface/NavigationSchool.h"
0022 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0023 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0024 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0025 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0026 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0027
0028 namespace {
0029
0030   Surface::RotationType rotation(const GlobalVector& zDir) {
0031     GlobalVector zAxis = zDir.unit();
0032     GlobalVector yAxis(zAxis.y(), -zAxis.x(), 0);
0033     GlobalVector xAxis = yAxis.cross(zAxis);
0034     return Surface::RotationType(xAxis, yAxis, zAxis);
0035   }
0036
0037 }  // namespace
0038
0039 class MeasurementTrackerTest : public edm::one::EDAnalyzer<> {
0040 public:
0041   explicit MeasurementTrackerTest(const edm::ParameterSet&);
0042   ~MeasurementTrackerTest() override = default;
0043
0044 private:
0045   virtual void analyze(const edm::Event&, const edm::EventSetup&) override;
0046
0047   const edm::ESGetToken<MeasurementTracker, CkfComponentsRecord> measTkToken;
0048   const edm::ESGetToken<NavigationSchool, NavigationSchoolRecord> navSchoolToken;
0049   const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mfToken;
0050   const edm::ESGetToken<Propagator, TrackingComponentsRecord> propToken;
0051
0052   std::string theMeasurementTrackerName;
0053   std::string theNavigationSchoolName;
0054 };
0055
0056 MeasurementTrackerTest::MeasurementTrackerTest(const edm::ParameterSet& iConfig)
0057     : measTkToken(esConsumes()),
0058       navSchoolToken(esConsumes()),
0059       mfToken(esConsumes()),
0060       propToken(esConsumes(edm::ESInputTag("", "PropagatorWithMaterial"))),
0061       theMeasurementTrackerName(iConfig.getParameter<std::string>("measurementTracker")),
0062       theNavigationSchoolName(iConfig.getParameter<std::string>("navigationSchool")) {}
0063
0064 void MeasurementTrackerTest::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0065   using namespace edm;
0066
0067   //get the measurementtracker
0068   const MeasurementTracker* measurementTracker = &iSetup.getData(measTkToken);
0069   const NavigationSchool* navSchool = &iSetup.getData(navSchoolToken);
0070
0071   auto const& geom = *(TrackerGeometry const*)(*measurementTracker).geomTracker();
0072   auto const& searchGeom = *(*measurementTracker).geometricSearchTracker();
0073   auto const& dus = geom.detUnits();
0074
0075   auto firstBarrel = geom.offsetDU(GeomDetEnumerators::tkDetEnum[1]);
0076   auto firstForward = geom.offsetDU(GeomDetEnumerators::tkDetEnum[2]);
0077
0078   std::cout << "number of dets " << dus.size() << std::endl;
0079   std::cout << "Bl/Fw loc " << firstBarrel << '/' << firstForward << std::endl;
0080
0081   const MagneticField* magfield = &iSetup.getData(mfToken);
0082   const Propagator& ANprop = iSetup.getData(propToken);
0083
0084   // error (very very small)
0085   ROOT::Math::SMatrixIdentity id;
0086   AlgebraicSymMatrix55 C(id);
0087   C *= 1.e-16;
0088
0089   CurvilinearTrajectoryError err(C);
0090
0091   //use negative sigma=-3.0 in order to use a more conservative definition of isInside() for Bounds classes.
0092   Chi2MeasurementEstimator estimator(30., -3.0, 0.5, 2.0, 0.5, 1.e12);  // same as defauts....
0093
0094   KFUpdator kfu;
0095   LocalError he(0.01 * 0.01, 0, 0.02 * 0.02);
0096
0097   for (float tl = 0.1f; tl < 3.0f; tl += 0.5f) {
0098     float p = 1.0f;
0099     float phi = 0.1415f;
0100     float tanlmd = tl;  // 0.1f;
0101     auto sinth2 = 1.f / (1.f + tanlmd * tanlmd);
0102     auto sinth = std::sqrt(sinth2);
0103     auto costh = tanlmd * sinth;
0104
0105     std::cout << tanlmd << ' ' << sinth << ' ' << costh << std::endl;
0106
0107     GlobalVector startingMomentum(p * std::sin(phi) * sinth, p * std::cos(phi) * sinth, p * costh);
0108     float z = 0.1f;
0109     GlobalPoint startingPosition(0, 0, z);
0110
0111     // make TSOS happy
0112     //Define starting plane
0113     PlaneBuilder pb;
0114     auto rot = rotation(startingMomentum);
0115     auto startingPlane = pb.plane(startingPosition, rot);
0116
0117     GlobalVector moms[3] = {0.5f * startingMomentum, startingMomentum, 10.f * startingMomentum};
0118
0119     for (auto mom : moms) {
0120       TrajectoryStateOnSurface startingStateP(
0121           GlobalTrajectoryParameters(startingPosition, mom, 1, magfield), err, *startingPlane);
0122       auto tsos = startingStateP;
0123
0124       // auto layer = searchGeom.idToLayer(dus[firstBarrel]->geographicalId());
0125       const DetLayer* layer = searchGeom.pixelBarrelLayers().front();
0126       {
0127         auto it = layer;
0128         std::cout << "first layer " << (it->isBarrel() ? " Barrel" : " Forward") << " layer " << it->seqNum()
0129                   << " SubDet " << it->subDetector() << std::endl;
0130       }
0131       auto const& detWithState = layer->compatibleDets(tsos, ANprop, estimator);
0132       if (!detWithState.size()) {
0133         std::cout << "no det on first layer" << std::endl;
0134         continue;
0135       }
0136       auto did = detWithState.front().first->geographicalId();
0137       std::cout << "arrived at " << int(did) << std::endl;
0138       tsos = detWithState.front().second;
0139       std::cout << tsos.globalPosition() << ' ' << tsos.localError().positionError() << std::endl;
0140
0141       SiPixelRecHit::ClusterRef pref;
0142       SiPixelRecHit hitpx(tsos.localPosition(), he, 1., *detWithState.front().first, pref);
0143       tsos = kfu.update(tsos, hitpx);
0144       std::cout << tsos.globalPosition() << ' ' << tsos.localError().positionError() << std::endl;
0145
0146       for (auto il = 1; il < 5; ++il) {
0147         if (!layer)
0148           break;
0149         auto const& compLayers = (*navSchool).nextLayers(*layer, *tsos.freeState(), alongMomentum);
0150         layer = nullptr;
0151         for (auto it : compLayers) {
0152           if (it->basicComponents().empty()) {
0153             //this should never happen. but better protect for it
0154             std::cout
0155                 << "a detlayer with no components: I cannot figure out a DetId from this layer. please investigate."
0156                 << std::endl;
0157             continue;
0158           }
0159           std::cout << il << (it->isBarrel() ? " Barrel" : " Forward") << " layer " << it->seqNum() << " SubDet "
0160                     << it->subDetector() << std::endl;
0161           auto const& detWithState = it->compatibleDets(tsos, ANprop, estimator);
0162           if (!detWithState.size()) {
0163             std::cout << "no det on this layer" << std::endl;
0164             continue;
0165           }
0166           layer = it;
0167           auto did = detWithState.front().first->geographicalId();
0168           std::cout << "arrived at " << int(did) << std::endl;
0169           tsos = detWithState.front().second;
0170           std::cout << tsos.globalPosition() << ' ' << tsos.localError().positionError() << std::endl;
0171         }
0172       }  // layer loop
0173     }  // loop on moms
0174   }  // loop  on tanLa
0175 }
0176
0177 //define this as a plug-in
0178 DEFINE_FWK_MODULE(MeasurementTrackerTest);