Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:29:00 

0001 #include <sstream>
0002 
0003 #include "DataFormats/Common/interface/OrphanHandle.h"
0004 #include "DataFormats/GeometryCommonDetAlgo/interface/ErrorFrameTransformer.h"
0005 #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
0006 #include "DataFormats/TrackReco/interface/TrackBase.h"
0007 #include "DataFormats/TrackerRecHit2D/interface/TrackingRecHitLessFromGlobalPosition.h"
0008 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0009 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012 #include "FWCore/Utilities/interface/thread_safety_macros.h"
0013 #include "Geometry/CommonDetUnit/interface/TrackingGeometry.h"
0014 #include "MagneticField/Engine/interface/MagneticField.h"
0015 #include "RecoTracker/TrackProducer/interface/TrackProducerAlgorithm.h"
0016 #include "RecoTracker/TransientTrackingRecHit/interface/TRecHit1DMomConstraint.h"
0017 #include "RecoTracker/TransientTrackingRecHit/interface/TRecHit2DPosConstraint.h"
0018 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0019 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
0020 #include "TrackingTools/PatternTools/interface/TSCBLBuilderWithPropagator.h"
0021 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0022 #include "TrackingTools/TrackFitters/interface/RecHitSorter.h"
0023 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0024 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0025 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0026 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0027 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHitBuilder.h"
0028 #include "TrackingTools/GsfTools/interface/GetComponents.h"
0029 
0030 // #define VI_DEBUG
0031 // #define STAT_TSB
0032 
0033 #ifdef VI_DEBUG
0034 #define DPRINT(x) std::cout << x << ": "
0035 #else
0036 #define DPRINT(x) LogTrace(x)
0037 #endif
0038 
0039 namespace {
0040 #ifdef STAT_TSB
0041   struct StatCount {
0042     long long totTrack = 0;
0043     long long totLoop = 0;
0044     long long totGsfTrack = 0;
0045     long long totFound = 0;
0046     long long totLost = 0;
0047     long long totAlgo[15];
0048     void track(int l) {
0049       if (l > 0)
0050         ++totLoop;
0051       else
0052         ++totTrack;
0053     }
0054     void hits(int f, int l) {
0055       totFound += f;
0056       totLost += l;
0057     }
0058     void gsf() { ++totGsfTrack; }
0059     void algo(int a) {
0060       if (a >= 0 && a < 15)
0061         ++totAlgo[a];
0062     }
0063 
0064     void print() const {
0065       std::cout << "TrackProducer stat\nTrack/Loop/Gsf/FoundHits/LostHits//algos " << totTrack << '/' << totLoop << '/'
0066                 << totGsfTrack << '/' << totFound << '/' << totLost << '/';
0067       for (auto a : totAlgo)
0068         std::cout << '/' << a;
0069       std::cout << std::endl;
0070     }
0071     StatCount() {}
0072     ~StatCount() { print(); }
0073   };
0074   StatCount statCount;
0075 
0076 #else
0077   struct StatCount {
0078     void track(int) {}
0079     void hits(int, int) {}
0080     void gsf() {}
0081     void algo(int) {}
0082   };
0083   CMS_THREAD_SAFE StatCount statCount;
0084 #endif
0085 
0086 }  // namespace
0087 
0088 template <>
0089 bool TrackProducerAlgorithm<reco::Track>::buildTrack(const TrajectoryFitter* theFitter,
0090                                                      const Propagator* thePropagator,
0091                                                      AlgoProductCollection& algoResults,
0092                                                      TransientTrackingRecHit::RecHitContainer& hits,
0093                                                      TrajectoryStateOnSurface& theTSOS,
0094                                                      const TrajectorySeed& seed,
0095                                                      float ndof,
0096                                                      const reco::BeamSpot& bs,
0097                                                      SeedRef seedRef,
0098                                                      int qualityMask,
0099                                                      signed char nLoops) {
0100   //variable declarations
0101 
0102   PropagationDirection seedDir = seed.direction();
0103 
0104   //perform the fit: the result's size is 1 if it succeded, 0 if fails
0105   Trajectory&& trajTmp =
0106       theFitter->fitOne(seed, hits, theTSOS, (nLoops > 0) ? TrajectoryFitter::looper : TrajectoryFitter::standard);
0107   if UNLIKELY (!trajTmp.isValid()) {
0108     DPRINT("TrackFitters") << "fit failed " << algo_ << ": " << hits.size() << '|' << int(nLoops) << ' ' << std::endl;
0109     return false;
0110   }
0111 
0112   auto theTraj = new Trajectory(std::move(trajTmp));
0113   theTraj->setSeedRef(seedRef);
0114 
0115   statCount.hits(theTraj->foundHits(), theTraj->lostHits());
0116   statCount.algo(int(algo_));
0117 
0118   // TrajectoryStateOnSurface innertsos;
0119   // if (theTraj->direction() == alongMomentum) {
0120   //  innertsos = theTraj->firstMeasurement().updatedState();
0121   // } else {
0122   //  innertsos = theTraj->lastMeasurement().updatedState();
0123   // }
0124 
0125   ndof = 0;
0126   for (auto const& tm : theTraj->measurements()) {
0127     auto const& h = tm.recHitR();
0128     if (h.isValid())
0129       ndof = ndof + float(h.dimension()) * h.weight();  // two virtual calls!
0130   }
0131 
0132   ndof -= 5.f;
0133   if UNLIKELY (theTSOS.magneticField()->nominalValue() == 0)
0134     ++ndof;  // same as -4
0135 
0136 #if defined(VI_DEBUG) || defined(EDM_ML_DEBUG)
0137   int chit[7] = {};
0138   int kk = 0;
0139   for (auto const& tm : theTraj->measurements()) {
0140     ++kk;
0141     auto const& hit = tm.recHitR();
0142     if (!hit.isValid())
0143       ++chit[0];
0144     if (hit.det() == nullptr)
0145       ++chit[1];
0146     if (trackerHitRTTI::isUndef(hit))
0147       continue;
0148     if (0)
0149       std::cout << "h " << kk << ": " << hit.localPosition() << ' ' << hit.localPositionError() << ' ' << tm.estimate()
0150                 << std::endl;
0151     if (hit.dimension() != 2) {
0152       ++chit[2];
0153     } else if (trackerHitRTTI::isFromDet(hit)) {
0154       auto const& thit = static_cast<BaseTrackerRecHit const&>(hit);
0155       auto const& clus = thit.firstClusterRef();
0156       if (clus.isPixel())
0157         ++chit[3];
0158       else if (thit.isMatched()) {
0159         ++chit[4];
0160       } else if (thit.isProjected()) {
0161         ++chit[5];
0162       } else {
0163         ++chit[6];
0164       }
0165     }
0166   }
0167 
0168   std::ostringstream ss;
0169   ss << algo_ << ": " << hits.size() << '|' << theTraj->measurements().size() << '|' << int(nLoops) << ' ';
0170   for (auto c : chit)
0171     ss << c << '/';
0172   ss << std::endl;
0173   DPRINT("TrackProducer") << ss.str();
0174 
0175 #endif
0176 
0177   //if geometricInnerState_ is false the state for projection to beam line is the state attached to the first hit: to be used for loopers
0178   //if geometricInnerState_ is true the state for projection to beam line is the one from the (geometrically) closest measurement to the beam line: to be sued for non-collision tracks
0179   //the two shouuld give the same result for collision tracks that are NOT loopers
0180   TrajectoryStateOnSurface stateForProjectionToBeamLineOnSurface;
0181   if (geometricInnerState_) {
0182     stateForProjectionToBeamLineOnSurface =
0183         theTraj->closestMeasurement(GlobalPoint(bs.x0(), bs.y0(), bs.z0())).updatedState();
0184   } else {
0185     if (theTraj->direction() == alongMomentum) {
0186       stateForProjectionToBeamLineOnSurface = theTraj->firstMeasurement().updatedState();
0187     } else {
0188       stateForProjectionToBeamLineOnSurface = theTraj->lastMeasurement().updatedState();
0189     }
0190   }
0191 
0192   if UNLIKELY (!stateForProjectionToBeamLineOnSurface.isValid()) {
0193     edm::LogError("CannotPropagateToBeamLine") << "the state on the closest measurement isnot valid. skipping track.";
0194     delete theTraj;
0195     return false;
0196   }
0197   const FreeTrajectoryState& stateForProjectionToBeamLine = *stateForProjectionToBeamLineOnSurface.freeState();
0198 
0199   LogDebug("TrackProducer") << "stateForProjectionToBeamLine=" << stateForProjectionToBeamLine;
0200 
0201   //  TSCBLBuilderNoMaterial tscblBuilder;
0202   //  TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,bs);
0203 
0204   TrajectoryStateClosestToBeamLine tscbl;
0205   if (usePropagatorForPCA_) {
0206     //std::cout << "PROPAGATOR FOR PCA" << std::endl;
0207     TSCBLBuilderWithPropagator tscblBuilder(*thePropagator);
0208     tscbl = tscblBuilder(stateForProjectionToBeamLine, bs);
0209   } else {
0210     TSCBLBuilderNoMaterial tscblBuilder;
0211     tscbl = tscblBuilder(stateForProjectionToBeamLine, bs);
0212   }
0213 
0214   if UNLIKELY (!tscbl.isValid()) {
0215     delete theTraj;
0216     return false;
0217   }
0218 
0219   GlobalPoint v = tscbl.trackStateAtPCA().position();
0220   math::XYZPoint pos(v.x(), v.y(), v.z());
0221   GlobalVector p = tscbl.trackStateAtPCA().momentum();
0222   math::XYZVector mom(p.x(), p.y(), p.z());
0223 
0224   LogDebug("TrackProducer") << "pos=" << v << " mom=" << p << " pt=" << p.perp() << " mag=" << p.mag();
0225 
0226   auto theTrack = new reco::Track(theTraj->chiSquared(),
0227                                   int(ndof),  //FIXME fix weight() in TrackingRecHit
0228                                   pos,
0229                                   mom,
0230                                   tscbl.trackStateAtPCA().charge(),
0231                                   tscbl.trackStateAtPCA().curvilinearError(),
0232                                   algo_);
0233 
0234   if (originalAlgo_ != reco::TrackBase::undefAlgorithm)
0235     theTrack->setOriginalAlgorithm(originalAlgo_);
0236   if (algoMask_.any())
0237     theTrack->setAlgoMask(algoMask_);
0238   theTrack->setQualityMask(qualityMask);
0239   theTrack->setNLoops(nLoops);
0240   theTrack->setStopReason(stopReason_);
0241 
0242   LogDebug("TrackProducer") << "theTrack->pt()=" << theTrack->pt();
0243 
0244   LogDebug("TrackProducer") << "track done\n";
0245 
0246   AlgoProduct aProduct{theTraj, theTrack, seedDir, 0};
0247   algoResults.push_back(aProduct);
0248 
0249   statCount.track(nLoops);
0250 
0251   return true;
0252 }
0253 
0254 template <>
0255 bool TrackProducerAlgorithm<reco::GsfTrack>::buildTrack(const TrajectoryFitter* theFitter,
0256                                                         const Propagator* thePropagator,
0257                                                         AlgoProductCollection& algoResults,
0258                                                         TransientTrackingRecHit::RecHitContainer& hits,
0259                                                         TrajectoryStateOnSurface& theTSOS,
0260                                                         const TrajectorySeed& seed,
0261                                                         float ndof,
0262                                                         const reco::BeamSpot& bs,
0263                                                         SeedRef seedRef,
0264                                                         int qualityMask,
0265                                                         signed char nLoops) {
0266   PropagationDirection seedDir = seed.direction();
0267 
0268   Trajectory&& trajTmp =
0269       theFitter->fitOne(seed, hits, theTSOS, (nLoops > 0) ? TrajectoryFitter::looper : TrajectoryFitter::standard);
0270   if UNLIKELY (!trajTmp.isValid())
0271     return false;
0272 
0273   auto theTraj = new Trajectory(std::move(trajTmp));
0274   theTraj->setSeedRef(seedRef);
0275 
0276 #ifdef EDM_ML_DEBUG
0277   TrajectoryStateOnSurface innertsos;
0278   TrajectoryStateOnSurface outertsos;
0279 
0280   if (theTraj->direction() == alongMomentum) {
0281     innertsos = theTraj->firstMeasurement().updatedState();
0282     outertsos = theTraj->lastMeasurement().updatedState();
0283   } else {
0284     innertsos = theTraj->lastMeasurement().updatedState();
0285     outertsos = theTraj->firstMeasurement().updatedState();
0286   }
0287   std::ostringstream ss;
0288   auto dc = [&](TrajectoryStateOnSurface const& tsos) {
0289     GetComponents comps(tsos);
0290     auto const& components = comps();
0291     auto sinTheta = std::sin(tsos.globalMomentum().theta());
0292     for (auto const& ic : components)
0293       ss << ic.weight() << "/";
0294     ss << "\n";
0295     for (auto const& ic : components)
0296       ss << ic.localParameters().vector()[0] / sinTheta << "/";
0297     ss << "\n";
0298     for (auto const& ic : components)
0299       ss << std::sqrt(ic.localError().matrix()(0, 0)) / sinTheta << "/";
0300   };
0301   ss << "\ninner comps\n";
0302   dc(innertsos);
0303   GetComponents icomps(innertsos);
0304   auto const& tsosComponentsInner = icomps();
0305 
0306   ss << "\nouter comps\n";
0307   dc(outertsos);
0308   GetComponents ocomps(outertsos);
0309   auto const& tsosComponentsOuter = ocomps();
0310 
0311   LogDebug("TrackProducer") << "Nr. of first / last states = " << tsosComponentsInner.size() << " "
0312                             << tsosComponentsOuter.size() << ss.str();
0313 #endif
0314 
0315   ndof = 0;
0316   for (auto const& tm : theTraj->measurements()) {
0317     auto const& h = tm.recHitR();
0318     if (h.isValid())
0319       ndof = ndof + h.dimension() * h.weight();
0320   }
0321 
0322   ndof = ndof - 5;
0323   if UNLIKELY (theTSOS.magneticField()->nominalValue() == 0)
0324     ++ndof;  // same as -4
0325 
0326   //if geometricInnerState_ is false the state for projection to beam line is the state attached to the first hit: to be used for loopers
0327   //if geometricInnerState_ is true the state for projection to beam line is the one from the (geometrically) closest measurement to the beam line: to be sued for non-collision tracks
0328   //the two shouuld give the same result for collision tracks that are NOT loopers
0329   TrajectoryStateOnSurface stateForProjectionToBeamLineOnSurface;
0330   if (geometricInnerState_) {
0331     stateForProjectionToBeamLineOnSurface =
0332         theTraj->closestMeasurement(GlobalPoint(bs.x0(), bs.y0(), bs.z0())).updatedState();
0333   } else {
0334     if (theTraj->direction() == alongMomentum) {
0335       stateForProjectionToBeamLineOnSurface = theTraj->firstMeasurement().updatedState();
0336     } else {
0337       stateForProjectionToBeamLineOnSurface = theTraj->lastMeasurement().updatedState();
0338     }
0339   }
0340 
0341   if UNLIKELY (!stateForProjectionToBeamLineOnSurface.isValid()) {
0342     edm::LogError("CannotPropagateToBeamLine") << "the state on the closest measurement isnot valid. skipping track.";
0343     delete theTraj;
0344     return false;
0345   }
0346 
0347   const FreeTrajectoryState& stateForProjectionToBeamLine = *stateForProjectionToBeamLineOnSurface.freeState();
0348 
0349   LogDebug("GsfTrackProducer") << "stateForProjectionToBeamLine=" << stateForProjectionToBeamLine;
0350 
0351   //  TSCBLBuilderNoMaterial tscblBuilder;
0352   //  TrajectoryStateClosestToBeamLine tscbl = tscblBuilder(stateForProjectionToBeamLine,bs);
0353 
0354   TrajectoryStateClosestToBeamLine tscbl;
0355   if (usePropagatorForPCA_) {
0356     TSCBLBuilderWithPropagator tscblBuilder(*thePropagator);
0357     tscbl = tscblBuilder(stateForProjectionToBeamLine, bs);
0358   } else {
0359     TSCBLBuilderNoMaterial tscblBuilder;
0360     tscbl = tscblBuilder(stateForProjectionToBeamLine, bs);
0361   }
0362 
0363   if UNLIKELY (tscbl.isValid() == false) {
0364     delete theTraj;
0365     return false;
0366   }
0367 
0368   GlobalPoint v = tscbl.trackStateAtPCA().position();
0369   math::XYZPoint pos(v.x(), v.y(), v.z());
0370   GlobalVector p = tscbl.trackStateAtPCA().momentum();
0371   math::XYZVector mom(p.x(), p.y(), p.z());
0372 
0373   LogDebug("GsfTrackProducer") << "pos=" << v << " mom=" << p << " pt=" << p.perp() << " mag=" << p.mag();
0374 
0375   auto theTrack =
0376       new reco::GsfTrack(theTraj->chiSquared(),
0377                          int(ndof),  //FIXME fix weight() in TrackingRecHit
0378                          //                theTraj->foundHits(),//FIXME to be fixed in Trajectory.h
0379                          //                0, //FIXME no corresponding method in trajectory.h
0380                          //                theTraj->lostHits(),//FIXME to be fixed in Trajectory.h
0381                          pos,
0382                          mom,
0383                          tscbl.trackStateAtPCA().charge(),
0384                          tscbl.trackStateAtPCA().curvilinearError());
0385   theTrack->setAlgorithm(algo_);
0386   if (originalAlgo_ != reco::TrackBase::undefAlgorithm)
0387     theTrack->setOriginalAlgorithm(originalAlgo_);
0388   if (algoMask_.any())
0389     theTrack->setAlgoMask(algoMask_);
0390 
0391   theTrack->setStopReason(stopReason_);
0392 
0393   LogDebug("GsfTrackProducer") << "track done\n";
0394 
0395   AlgoProduct aProduct{theTraj, theTrack, seedDir, 0};
0396 
0397   LogDebug("GsfTrackProducer") << "track done1\n";
0398   algoResults.push_back(aProduct);
0399   LogDebug("GsfTrackProducer") << "track done2\n";
0400 
0401   statCount.gsf();
0402   return true;
0403 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!