TkNavigation/plugins/TkMSParameterizationBuilder.cc

0001 #include "RecoTracker/TkNavigation/interface/TkMSParameterization.h"
0002
0003 #include "RecoTracker/TkNavigation/interface/TkNavigationSchool.h"
0004 #include "RecoTracker/Record/interface/NavigationSchoolRecord.h"
0005 #include "MagneticField/Engine/interface/MagneticField.h"
0006
0007 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0008 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0009 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0010 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0011
0012 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0013 #include "TrackPropagation/RungeKutta/interface/defaultRKPropagator.h"
0014 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0015
0016 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0017 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0018
0019 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0020
0021 #include "FWCore/Framework/interface/ESProducer.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0025
0026 #include "RecoTracker/Record/interface/TkMSParameterizationRecord.h"
0027
0028 namespace {
0029   inline Surface::RotationType rotation(const GlobalVector& zDir) {
0030     GlobalVector zAxis = zDir.unit();
0031     GlobalVector yAxis(zAxis.y(), -zAxis.x(), 0);
0032     GlobalVector xAxis = yAxis.cross(zAxis);
0033     return Surface::RotationType(xAxis, yAxis, zAxis);
0034   }
0035
0036   struct MSData {
0037     int stid;
0038     int lid;
0039     float zi;
0040     float zo;
0041     float uerr;
0042     float verr;
0043   };
0044 }  // namespace
0045 inline std::ostream& operator<<(std::ostream& os, MSData d) {
0046   os << d.stid << '>' << d.lid << '|' << d.zi << '/' << d.zo << ':' << d.uerr << '/' << d.verr;
0047   return os;
0048 }
0049
0050 class TkMSParameterizationBuilder final : public edm::ESProducer {
0051 public:
0052   TkMSParameterizationBuilder(edm::ParameterSet const&);
0053   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0054     edm::ParameterSetDescription desc;
0055     // desc.add<std::string>("ComponentName","TkMSParameterization");
0056     desc.add<std::string>("navigationSchool", "SimpleNavigationSchool");
0057     descriptions.add("TkMSParameterizationBuilder", desc);
0058   }
0059
0060   using ReturnType = std::unique_ptr<TkMSParameterization>;
0061   ReturnType produce(TkMSParameterizationRecord const&);
0062
0063 private:
0064   edm::ESGetToken<NavigationSchool, NavigationSchoolRecord> theNavSchoolToken_;
0065   edm::ESGetToken<Propagator, TrackingComponentsRecord> thePropagatorToken_;
0066 };
0067
0068 TkMSParameterizationBuilder::TkMSParameterizationBuilder(edm::ParameterSet const& pset) {
0069   auto cc = setWhatProduced(this, "");
0070   theNavSchoolToken_ = cc.consumes(edm::ESInputTag("", pset.getParameter<std::string>("navigationSchool")));
0071   thePropagatorToken_ = cc.consumes(edm::ESInputTag("", "PropagatorWithMaterial"));
0072 }
0073
0074 TkMSParameterizationBuilder::ReturnType TkMSParameterizationBuilder::produce(TkMSParameterizationRecord const& iRecord) {
0075   using namespace tkMSParameterization;
0076
0077   ReturnType product = std::make_unique<TkMSParameterization>();
0078
0079   auto& msParam = *product;
0080
0081   //
0082   TkNavigationSchool const& navSchool = dynamic_cast<TkNavigationSchool const&>(iRecord.get(theNavSchoolToken_));
0083   auto const& searchGeom = navSchool.searchTracker();
0084   auto const& magfield = navSchool.field();
0085
0086   //
0087   auto const& ANprop = iRecord.get(thePropagatorToken_);
0088
0089   // error (very very small)
0090   ROOT::Math::SMatrixIdentity id;
0091   AlgebraicSymMatrix55 C(id);
0092   C *= 1.e-16;
0093
0094   CurvilinearTrajectoryError err(C);
0095
0096   Chi2MeasurementEstimator estimator(30., -3.0, 0.5, 2.0, 0.5, 1.e12);  // same as defauts....
0097   KFUpdator kfu;
0098   LocalError he(0.001 * 0.001, 0, 0.002 * 0.002);
0099
0100   // loop over lambdas
0101   bool debug = false;
0102   float tl = 0;
0103   for (decltype(nLmBins()) ib = 0; ib < nLmBins(); ++ib, tl += lmBin()) {
0104     float pt = 1.0f;
0105     float tanlmd = tl;  // 0.1f;
0106     auto sinth2 = 1.f / (1.f + tanlmd * tanlmd);
0107     auto sinth = std::sqrt(sinth2);
0108     auto overp = sinth / pt;
0109     auto pz = pt * tanlmd;
0110
0111     // debug= (tl>2.34f && tl<2.55f);
0112     if (debug)
0113       std::cout << tl << ' ' << pz << ' ' << 1. / overp << std::endl;
0114
0115     float lastzz = -18.f;
0116     bool goFw = false;
0117     std::string loc = " Barrel";
0118     for (int iz = 0; iz < 2; ++iz) {
0119       if (iz > 0)
0120         goFw = true;
0121       for (float zz = lastzz; zz < 18.1f; zz += 0.2f) {
0122         float z = zz;
0123         GlobalPoint startingPosition(0, 0, z);
0124
0125         constexpr int maxLayers = 6;
0126
0127         std::vector<MSData> mserr[maxLayers][maxLayers];
0128
0129         // define propagation from/to
0130         std::function<void(int, TrajectoryStateOnSurface, DetLayer const*, float, int)> propagate =
0131             [&](int from, TrajectoryStateOnSurface tsos, DetLayer const* layer, float z1, int stid) {
0132               for (auto il = from + 1; il < maxLayers; ++il) {
0133                 auto compLayers = navSchool.nextLayers(*layer, *tsos.freeState(), alongMomentum);
0134                 std::stable_sort(
0135                     compLayers.begin(), compLayers.end(), [](auto a, auto b) { return a->seqNum() < b->seqNum(); });
0136                 layer = nullptr;
0137                 for (auto it : compLayers) {
0138                   if (it->basicComponents().empty()) {
0139                     //this should never happen. but better protect for it
0140                     edm::LogError("TkMSParameterizationBuilder") << "a detlayer with no components: I cannot figure "
0141                                                                     "out a DetId from this layer. please investigate.";
0142                     continue;
0143                   }
0144                   if (debug)
0145                     std::cout << il << (it->isBarrel() ? " Barrel" : " Forward") << " layer " << it->seqNum()
0146                               << " SubDet " << it->subDetector() << std::endl;
0147                   auto const& detWithState = it->compatibleDets(tsos, ANprop, estimator);
0148                   if (detWithState.empty()) {
0149                     if (debug)
0150                       std::cout << "no det on this layer" << it->seqNum() << std::endl;
0151                     continue;
0152                   }
0153                   layer = it;
0154                   auto did = detWithState.front().first->geographicalId();
0155                   if (debug)
0156                     std::cout << "arrived at " << int(did) << std::endl;
0157                   tsos = detWithState.front().second;
0158                   if (debug)
0159                     std::cout << tsos.globalPosition() << ' ' << tsos.localError().positionError() << std::endl;
0160
0161                   // save errs
0162                   auto globalError =
0163                       ErrorFrameTransformer::transform(tsos.localError().positionError(), tsos.surface());
0164                   auto gp = tsos.globalPosition();
0165                   float r = gp.perp();
0166                   float errorPhi = std::sqrt(float(globalError.phierr(gp)));
0167                   float errorR = std::sqrt(float(globalError.rerr(gp)));
0168                   float errorZ = std::sqrt(float(globalError.czz()));
0169                   float xerr = overp * errorPhi;
0170                   float zerr = overp * (layer->isBarrel() ? errorZ : errorR);
0171                   auto zo = layer->isBarrel() ? gp.z() : r;
0172                   //  std::cout << tanlmd << ' ' << z1 << ' ' << it->seqNum() << ':' << xerr <<'/'<<zerr << std::endl;
0173                   mserr[from][il - 1].emplace_back(MSData{stid, it->seqNum(), z1, zo, xerr, zerr});
0174
0175                   if (from == 0) {
0176                     // spawn prop from this layer
0177                     // constrain it to this location (relevant for layer other than very first)
0178                     SiPixelRecHit::ClusterRef pref;
0179                     SiPixelRecHit hitpx(tsos.localPosition(), he, 1., *detWithState.front().first, pref);
0180                     auto tsosl = kfu.update(tsos, hitpx);
0181                     auto z1l = layer->isBarrel() ? tsos.globalPosition().z() : tsos.globalPosition().perp();
0182                     auto stidl = layer->seqNum();
0183                     propagate(il, tsosl, layer, z1l, stidl);
0184                   }
0185                   break;
0186                 }  // loop on compLayers
0187                 if (!layer)
0188                   break;
0189                 // if (!goFw) lastbz=z1;
0190                 if (from == 0)
0191                   lastzz = zz;
0192
0193               }  // layer loop
0194
0195               if (debug && mserr[from][from].empty()) {
0196                 std::cout << "tl " << tanlmd << loc << ' ' << from << std::endl;
0197                 for (auto il = from; il < maxLayers; ++il) {
0198                   std::cout << il << ' ';
0199                   for (auto const& e : mserr[from][il])
0200                     std::cout << e << ' ';  //  << '-' <<e.uerr*sinth <<'/'<<e.verr*sinth <<' ';
0201                   std::cout << std::endl;
0202                 }
0203               }
0204             };  // propagate
0205
0206         float phi = 0.0415f;
0207         for (int ip = 0; ip < 5; ++ip) {
0208           phi += (3.14159f / 6.f) / 5.f;
0209           GlobalVector startingMomentum(pt * std::sin(phi), pt * std::cos(phi), pz);
0210
0211           // make TSOS happy
0212           //Define starting plane
0213           PlaneBuilder pb;
0214           auto rot = rotation(startingMomentum);
0215           auto startingPlane = pb.plane(startingPosition, rot);
0216
0217           TrajectoryStateOnSurface startingStateP(
0218               GlobalTrajectoryParameters(startingPosition, startingMomentum, 1, &magfield), err, *startingPlane);
0219           auto tsos0 = startingStateP;
0220
0221           DetLayer const* layer0 = searchGeom.pixelBarrelLayers()[0];
0222           if (goFw)
0223             layer0 = searchGeom.posPixelForwardLayers()[0];
0224           int stid0 = layer0->seqNum();
0225
0226           {
0227             auto it = layer0;
0228             if (debug)
0229               std::cout << "first layer " << (it->isBarrel() ? " Barrel" : " Forward") << " layer " << it->seqNum()
0230                         << " SubDet " << it->subDetector() << std::endl;
0231           }
0232
0233           auto const& detWithState = layer0->compatibleDets(tsos0, ANprop, estimator);
0234           if (detWithState.empty()) {
0235             if (debug)
0236               std::cout << "no det on first layer" << layer0->seqNum() << std::endl;
0237             continue;
0238           }
0239           tsos0 = detWithState.front().second;
0240           if (debug)
0241             std::cout << "arrived at " << int(detWithState.front().first->geographicalId()) << ' '
0242                       << tsos0.globalPosition() << ' ' << tsos0.localError().positionError() << std::endl;
0243
0244           // for barrel
0245           float z1l = tsos0.globalPosition().z();
0246           if (goFw) {
0247             loc = " Forward";
0248             z1l = tsos0.globalPosition().perp();
0249           }
0250
0251           propagate(0, tsos0, layer0, z1l, stid0);
0252
0253         }  // loop on phi
0254
0255         // fill
0256         for (auto from = 0; from < maxLayers; ++from)
0257           for (auto il = from; il < maxLayers; ++il) {
0258             if (mserr[from][il].empty())
0259               continue;
0260             auto stid = mserr[from][il].front().stid;
0261             auto lid = mserr[from][il].front().lid;
0262             auto zi = mserr[from][il].front().zi;
0263             float xerr = 0;
0264             float zerr = 0;
0265             float zo = 0;
0266             for (auto const& e : mserr[from][il]) {
0267               if (e.stid != stid)
0268                 continue;
0269               lid = std::min(lid, e.lid);
0270             }
0271             float nn = 0;
0272             for (auto const& e : mserr[from][il]) {
0273               if (e.stid != stid || lid != e.lid)
0274                 continue;
0275               xerr += e.uerr;
0276               zerr += e.verr;
0277               zo += e.zo;
0278               ++nn;
0279             }
0280             xerr /= nn;
0281             zerr /= nn;
0282             zo /= nn;
0283             auto iid = packLID(stid, lid);
0284             auto& md = msParam.data[iid].data[ib].data;
0285             if (md.empty() || std::abs(xerr - md.back().uerr) > 0.2f * xerr ||
0286                 std::abs(zerr - md.back().verr) > 0.2f * zerr)
0287               md.emplace_back(Elem{zi, zo, xerr, zerr});
0288           }  // loop on layers
0289       }
0290     }  // loop on z
0291
0292   }  // loop  on tanLa
0293
0294   // sort
0295   for (auto& e : msParam.data)
0296     for (auto& d : e.second.data)
0297       std::stable_sort(d.data.begin(), d.data.end(), [](auto const& a, auto const& b) { return a.vo < b.vo; });
0298
0299   return product;
0300 }
0301
0302 #include "FWCore/PluginManager/interface/ModuleDef.h"
0303 #include "FWCore/Framework/interface/MakerMacros.h"
0304 #include "FWCore/Framework/interface/ModuleFactory.h"
0305 #include "FWCore/Utilities/interface/typelookup.h"
0306 DEFINE_FWK_EVENTSETUP_MODULE(TkMSParameterizationBuilder);