TkSeedGenerator/plugins/SeedingOTEDProducer.cc

0001
0002 #include "FWCore/Framework/interface/stream/EDProducer.h"
0003 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005
0006 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/VectorHit.h"
0008 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0009 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0010 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0011 #include "TrackingTools/MeasurementDet/interface/LayerMeasurements.h"
0012 #include "MagneticField/Engine/interface/MagneticField.h"
0013
0014 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0015 #include "RecoTracker/MeasurementDet/interface/MeasurementTracker.h"
0016 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0017 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0018 #include "TrackingTools/PatternTools/interface/TrajectoryStateUpdator.h"
0019
0020 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0021
0022 #include "FWCore/Utilities/interface/ESGetToken.h"
0023
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/MakerMacros.h"
0026 #include "FWCore/Framework/interface/EventSetup.h"
0027
0028 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0029 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0030 #include "TrackingTools/MeasurementDet/interface/TrajectoryMeasurementGroup.h"
0031
0032 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034
0035 class TrajectoryStateUpdator;
0036
0037 class SeedingOTEDProducer : public edm::stream::EDProducer<> {
0038 public:
0039   explicit SeedingOTEDProducer(const edm::ParameterSet&);
0040   ~SeedingOTEDProducer() override;
0041   void produce(edm::Event&, const edm::EventSetup&) override;
0042
0043   static void fillDescriptions(edm::ConfigurationDescriptions&);
0044
0045   TrajectorySeedCollection run(edm::Handle<VectorHitCollection>);
0046   unsigned int checkLayer(unsigned int iidd);
0047   std::vector<const VectorHit*> collectVHsOnLayer(const edmNew::DetSetVector<VectorHit>&, unsigned int);
0048   void printVHsOnLayer(const edmNew::DetSetVector<VectorHit>&, unsigned int);
0049   const TrajectoryStateOnSurface buildInitialTSOS(const VectorHit*) const;
0050   AlgebraicSymMatrix55 assign44To55(const AlgebraicSymMatrix44&) const;
0051   std::pair<bool, TrajectoryStateOnSurface> propagateAndUpdate(const TrajectoryStateOnSurface initialTSOS,
0052                                                                const Propagator&,
0053                                                                const TrackingRecHit& hit) const;
0054   float computeGlobalThetaError(const VectorHit* vh, const double sigmaZ_beamSpot) const;
0055   float computeInverseMomentumError(const VectorHit* vh,
0056                                     const float globalTheta,
0057                                     const double sigmaZ_beamSpot,
0058                                     const double transverseMomentum) const;
0059
0060   TrajectorySeed createSeed(const TrajectoryStateOnSurface& tsos,
0061                             const edm::OwnVector<TrackingRecHit>& container,
0062                             const DetId& id,
0063                             const Propagator& prop) const;
0064
0065   struct isInvalid {
0066     bool operator()(const TrajectoryMeasurement& measurement) {
0067       return ((measurement.recHit() == nullptr) || !(measurement.recHit()->isValid()) ||
0068               !((measurement).updatedState().isValid()));
0069     }
0070   };
0071
0072 private:
0073   edm::EDGetTokenT<VectorHitCollection> vhProducerToken_;
0074   const TrackerTopology* tkTopo_;
0075   const MeasurementTracker* measurementTracker_;
0076   std::unique_ptr<LayerMeasurements> layerMeasurements_;
0077   const MeasurementEstimator* estimator_;
0078   const Propagator* propagator_;
0079   const MagneticField* magField_;
0080   const TrajectoryStateUpdator* updator_;
0081   const edm::EDGetTokenT<MeasurementTrackerEvent> tkMeasEventToken_;
0082   edm::EDGetTokenT<reco::BeamSpot> beamSpotToken_;
0083   const reco::BeamSpot* beamSpot_;
0084   std::string updatorName_;
0085
0086   edm::ESGetToken<TrackerTopology, TrackerTopologyRcd> topoToken_;
0087   edm::ESGetToken<Propagator, TrackingComponentsRecord> propagatorToken_;
0088   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0089   edm::ESGetToken<TrajectoryStateUpdator, TrackingComponentsRecord> updatorToken_;
0090   edm::ESGetToken<MeasurementTracker, CkfComponentsRecord> measurementTrackerToken_;
0091   edm::ESGetToken<Chi2MeasurementEstimatorBase, TrackingComponentsRecord> estToken_;
0092   edm::EDPutTokenT<TrajectorySeedCollection> putToken_;
0093 };
0094
0095 SeedingOTEDProducer::SeedingOTEDProducer(edm::ParameterSet const& conf)
0096     : tkMeasEventToken_(consumes<MeasurementTrackerEvent>(conf.getParameter<edm::InputTag>("trackerEvent"))),
0097       topoToken_(esConsumes()),
0098       propagatorToken_(esConsumes(edm::ESInputTag("", "PropagatorWithMaterial"))),
0099       magFieldToken_(esConsumes()),
0100       updatorToken_(esConsumes(edm::ESInputTag("", "KFUpdator"))),
0101       measurementTrackerToken_(esConsumes()),
0102       estToken_(esConsumes(edm::ESInputTag("", "Chi2"))) {
0103   vhProducerToken_ = consumes<VectorHitCollection>(edm::InputTag(conf.getParameter<edm::InputTag>("src")));
0104   beamSpotToken_ = consumes<reco::BeamSpot>(conf.getParameter<edm::InputTag>("beamSpotLabel"));
0105   updatorName_ = conf.getParameter<std::string>("updator");
0106   putToken_ = produces<TrajectorySeedCollection>();
0107 }
0108
0109 SeedingOTEDProducer::~SeedingOTEDProducer() {}
0110
0111 void SeedingOTEDProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0112   edm::ParameterSetDescription desc;
0113   desc.add<edm::InputTag>("src", edm::InputTag("siPhase2VectorHits", "accepted"));
0114   desc.add<edm::InputTag>("trackerEvent", edm::InputTag("MeasurementTrackerEvent"));
0115   desc.add<edm::InputTag>("beamSpotLabel", edm::InputTag("offlineBeamSpot"));
0116   desc.add<std::string>("updator", std::string("KFUpdator"));
0117   descriptions.add("SeedingOTEDProducer", desc);
0118 }
0119
0120 void SeedingOTEDProducer::produce(edm::Event& event, const edm::EventSetup& es) {
0121   tkTopo_ = &es.getData(topoToken_);
0122
0123   edm::ESHandle<MeasurementTracker> measurementTrackerHandle = es.getHandle(measurementTrackerToken_);
0124   measurementTracker_ = measurementTrackerHandle.product();
0125
0126   auto const& measurementTrackerEvent = event.get(tkMeasEventToken_);
0127
0128   layerMeasurements_ = std::make_unique<LayerMeasurements>(*measurementTrackerHandle, measurementTrackerEvent);
0129
0130   estimator_ = &es.getData(estToken_);
0131
0132   propagator_ = &es.getData(propagatorToken_);
0133
0134   magField_ = &es.getData(magFieldToken_);
0135
0136   updator_ = &es.getData(updatorToken_);
0137
0138   beamSpot_ = &event.get(beamSpotToken_);
0139
0140   // Get the vector hits
0141   auto vhs = event.getHandle(vhProducerToken_);
0142
0143   event.emplace(putToken_, run(vhs));
0144 }
0145
0146 TrajectorySeedCollection SeedingOTEDProducer::run(edm::Handle<VectorHitCollection> vHs) {
0147   TrajectorySeedCollection result;
0148   //check if all the first three layers have VHs
0149   std::vector<const VectorHit*> vhSeedsL1 = collectVHsOnLayer(*vHs.product(), 1);
0150   std::vector<const VectorHit*> vhSeedsL2 = collectVHsOnLayer(*vHs.product(), 2);
0151   std::vector<const VectorHit*> vhSeedsL3 = collectVHsOnLayer(*vHs.product(), 3);
0152   if (vhSeedsL1.empty() || vhSeedsL2.empty() || vhSeedsL3.empty()) {
0153     return result;
0154   }
0155
0156   //seeds are built in the L3 of the OT
0157   const BarrelDetLayer* barrelOTLayer2 = measurementTracker_->geometricSearchTracker()->tobLayers().at(1);
0158
0159   //the search propag directiondepend on the sign of signZ*signPz, while the building is always the contrary
0160   Propagator* searchingPropagator = &*propagator_->clone();
0161   Propagator* buildingPropagator = &*propagator_->clone();
0162   buildingPropagator->setPropagationDirection(alongMomentum);
0163
0164   for (const auto& hitL3 : vhSeedsL3) {
0165     //building a tsos out of a VectorHit
0166     const TrajectoryStateOnSurface initialTSOS = buildInitialTSOS(hitL3);
0167     float signZ = copysign(1.0, initialTSOS.globalPosition().z());
0168     float signPz = copysign(1.0, initialTSOS.globalMomentum().z());
0169
0170     //set the direction of the propagator
0171     if (signZ * signPz > 0.0)
0172       searchingPropagator->setPropagationDirection(oppositeToMomentum);
0173     else
0174       searchingPropagator->setPropagationDirection(alongMomentum);
0175
0176     //find vHits in layer 2
0177     std::vector<TrajectoryMeasurement> measurementsL2 =
0178         layerMeasurements_->measurements(*barrelOTLayer2, initialTSOS, *searchingPropagator, *estimator_);
0179
0180     std::vector<TrajectoryMeasurement>::iterator measurementsL2end =
0181         std::remove_if(measurementsL2.begin(), measurementsL2.end(), isInvalid());
0182     measurementsL2.erase(measurementsL2end, measurementsL2.end());
0183
0184     if (!measurementsL2.empty()) {
0185       //not sure if building it everytime takes time/memory
0186       const DetLayer* barrelOTLayer1 = measurementTracker_->geometricSearchTracker()->tobLayers().at(0);
0187
0188       for (const auto& mL2 : measurementsL2) {
0189         const TrackingRecHit* hitL2 = mL2.recHit().get();
0190
0191         //propagate to the L2 and update the TSOS
0192         std::pair<bool, TrajectoryStateOnSurface> updatedTSOS =
0193             propagateAndUpdate(initialTSOS, *searchingPropagator, *hitL2);
0194         if (!updatedTSOS.first)
0195           continue;
0196
0197         //searching possible VHs in L1
0198         std::vector<TrajectoryMeasurement> measurementsL1 =
0199             layerMeasurements_->measurements(*barrelOTLayer1, updatedTSOS.second, *searchingPropagator, *estimator_);
0200         std::vector<TrajectoryMeasurement>::iterator measurementsL1end =
0201             std::remove_if(measurementsL1.begin(), measurementsL1.end(), isInvalid());
0202         measurementsL1.erase(measurementsL1end, measurementsL1.end());
0203
0204         if (!measurementsL1.empty()) {
0205           for (const auto& mL1 : measurementsL1) {
0206             const TrackingRecHit* hitL1 = mL1.recHit().get();
0207
0208             //propagate to the L1 and update the TSOS
0209             std::pair<bool, TrajectoryStateOnSurface> updatedTSOSL1 =
0210                 propagateAndUpdate(updatedTSOS.second, *searchingPropagator, *hitL1);
0211             if (!updatedTSOSL1.first)
0212               continue;
0213
0214             //building trajectory inside-out
0215             if (searchingPropagator->propagationDirection() == alongMomentum) {
0216               buildingPropagator->setPropagationDirection(oppositeToMomentum);
0217             } else {
0218               buildingPropagator->setPropagationDirection(alongMomentum);
0219             }
0220
0221             updatedTSOSL1.second.rescaleError(100);
0222
0223             TrajectoryStateOnSurface updatedTSOSL1_final = updator_->update(updatedTSOSL1.second, *hitL1);
0224             if UNLIKELY (!updatedTSOSL1_final.isValid())
0225               continue;
0226             std::pair<bool, TrajectoryStateOnSurface> updatedTSOSL2_final =
0227                 propagateAndUpdate(updatedTSOSL1_final, *buildingPropagator, *hitL2);
0228             if (!updatedTSOSL2_final.first)
0229               continue;
0230             std::pair<bool, TrajectoryStateOnSurface> updatedTSOSL3_final =
0231                 propagateAndUpdate(updatedTSOSL2_final.second, *buildingPropagator, *hitL3);
0232             if (!updatedTSOSL3_final.first)
0233               continue;
0234
0235             edm::OwnVector<TrackingRecHit> container;
0236             container.push_back(hitL1->clone());
0237             container.push_back(hitL2->clone());
0238             container.push_back(hitL3->clone());
0239
0240             TrajectorySeed ts =
0241                 createSeed(updatedTSOSL3_final.second, container, hitL3->geographicalId(), *buildingPropagator);
0242             result.push_back(ts);
0243           }
0244         }
0245       }
0246     }
0247   }
0248   result.shrink_to_fit();
0249   return result;
0250 }
0251
0252 unsigned int SeedingOTEDProducer::checkLayer(unsigned int iidd) {
0253   StripSubdetector strip = StripSubdetector(iidd);
0254   unsigned int subid = strip.subdetId();
0255   if (subid == StripSubdetector::TIB || subid == StripSubdetector::TOB) {
0256     return tkTopo_->layer(iidd);
0257   }
0258   return 0;
0259 }
0260
0261 std::vector<const VectorHit*> SeedingOTEDProducer::collectVHsOnLayer(const edmNew::DetSetVector<VectorHit>& input,
0262                                                                      unsigned int layerNumber) {
0263   std::vector<const VectorHit*> vHsOnLayer;
0264   if (!input.empty()) {
0265     for (const auto& DSViter : input) {
0266       if (checkLayer(DSViter.id()) == layerNumber) {
0267         for (const auto& vh : DSViter) {
0268           vHsOnLayer.push_back(&vh);
0269         }
0270       }
0271     }
0272   }
0273
0274   return vHsOnLayer;
0275 }
0276
0277 void SeedingOTEDProducer::printVHsOnLayer(const edmNew::DetSetVector<VectorHit>& input, unsigned int layerNumber) {
0278   if (!input.empty()) {
0279     for (const auto& DSViter : input) {
0280       for (const auto& vh : DSViter) {
0281         if (checkLayer(DSViter.id()) == layerNumber)
0282           LogTrace("SeedingOTEDProducer") << " VH in layer " << layerNumber << " >> " << vh;
0283       }
0284     }
0285   } else {
0286     LogTrace("SeedingOTEDProducer") << " No VHs in layer " << layerNumber << ".";
0287   }
0288 }
0289
0290 const TrajectoryStateOnSurface SeedingOTEDProducer::buildInitialTSOS(const VectorHit* vHit) const {
0291   // having fun with theta
0292   Global3DVector gv(vHit->globalPosition().x(), vHit->globalPosition().y(), vHit->globalPosition().z());
0293   float theta = gv.theta();
0294   // gv transform to local (lv)
0295   const Local3DVector lv(vHit->det()->surface().toLocal(gv));
0296
0297   //FIXME::charge is fine 1 every two times!!
0298   GlobalPoint center(0.0, 0.0, 0.0);
0299   int charge = 1;
0300   // momentum is a signed quantity in this case
0301   float mom = vHit->momentum(magField_->inTesla(center).z());
0302   float xPos = vHit->localPosition().x();
0303   float yPos = vHit->localPosition().y();
0304   float dx = vHit->localDirection().x();
0305   // for dy use second component of the lv renormalized to the z component
0306   float dy = lv.y() / lv.z();
0307
0308   // Pz and Dz should have the same sign
0309   float signPz = copysign(1.0, vHit->globalPosition().z());
0310
0311   LocalTrajectoryParameters ltpar2(charge / mom, dx, dy, xPos, yPos, signPz);
0312   AlgebraicSymMatrix55 mat = assign44To55(vHit->covMatrix());
0313   // set the error on 1/p
0314   mat[0][0] = pow(computeInverseMomentumError(
0315                       vHit, theta, beamSpot_->sigmaZ(), vHit->transverseMomentum(magField_->inTesla(center).z())),
0316                   2);
0317
0318   //building tsos
0319   LocalTrajectoryError lterr(mat);
0320   const TrajectoryStateOnSurface tsos(ltpar2, lterr, vHit->det()->surface(), magField_);
0321
0322   return tsos;
0323 }
0324
0325 AlgebraicSymMatrix55 SeedingOTEDProducer::assign44To55(const AlgebraicSymMatrix44& mat44) const {
0326   AlgebraicSymMatrix55 result;
0327   for (int i = 1; i < 5; i++) {
0328     for (int j = 1; j < 5; j++) {
0329       result[i][j] = mat44[i - 1][j - 1];
0330     }
0331   }
0332   return result;
0333 }
0334
0335 std::pair<bool, TrajectoryStateOnSurface> SeedingOTEDProducer::propagateAndUpdate(
0336     const TrajectoryStateOnSurface initialTSOS, const Propagator& prop, const TrackingRecHit& hit) const {
0337   TrajectoryStateOnSurface propTSOS = prop.propagate(initialTSOS, hit.det()->surface());
0338   if UNLIKELY (!propTSOS.isValid())
0339     return std::make_pair(false, propTSOS);
0340   TrajectoryStateOnSurface updatedTSOS = updator_->update(propTSOS, hit);
0341   if UNLIKELY (!updatedTSOS.isValid())
0342     return std::make_pair(false, updatedTSOS);
0343   return std::make_pair(true, updatedTSOS);
0344 }
0345
0346 float SeedingOTEDProducer::computeGlobalThetaError(const VectorHit* vh, const double sigmaZ_beamSpot) const {
0347   double derivative = vh->globalPosition().perp() / vh->globalPosition().mag2();
0348   double derivative2 = pow(derivative, 2);
0349   return pow(derivative2 * vh->lowerGlobalPosErr().czz() + derivative2 * pow(sigmaZ_beamSpot, 2), 0.5);
0350 }
0351
0352 float SeedingOTEDProducer::computeInverseMomentumError(const VectorHit* vh,
0353                                                        const float globalTheta,
0354                                                        const double sigmaZ_beamSpot,
0355                                                        const double transverseMomentum) const {
0356   //for pT > 2GeV, 1/pT has sigma = 1/sqrt(12)
0357   float varianceInverseTransvMomentum = 1. / 12;
0358   float derivativeTheta2 = pow(cos(globalTheta) / transverseMomentum, 2);
0359   float derivativeInverseTransvMomentum2 = pow(sin(globalTheta), 2);
0360   float thetaError = computeGlobalThetaError(vh, sigmaZ_beamSpot);
0361   return pow(derivativeTheta2 * pow(thetaError, 2) + derivativeInverseTransvMomentum2 * varianceInverseTransvMomentum,
0362              0.5);
0363 }
0364
0365 TrajectorySeed SeedingOTEDProducer::createSeed(const TrajectoryStateOnSurface& tsos,
0366                                                const edm::OwnVector<TrackingRecHit>& container,
0367                                                const DetId& id,
0368                                                const Propagator& prop) const {
0369   PTrajectoryStateOnDet seedTSOS = trajectoryStateTransform::persistentState(tsos, id.rawId());
0370   return TrajectorySeed(seedTSOS, container, prop.propagationDirection());
0371 }
0372 DEFINE_FWK_MODULE(SeedingOTEDProducer);