CMS 3D CMS Logo

Trajectory.cc
Go to the documentation of this file.
6 #include "boost/intrusive_ptr.hpp"
10 
11 
12 #include <algorithm>
13 
14 namespace {
15  template<typename DataContainer>
16  unsigned short countTrailingValidHits(DataContainer const & meas) {
17  unsigned short n=0;
18  for(auto it=meas.rbegin(); it!=meas.rend(); ++it) {
19  if (Trajectory::lost(*(*it).recHit())) break;
20  if((*it).recHit()->isValid()) ++n;
21  }
22  return n;
23  }
24 
25 }
26 
27 using namespace std;
28 
30  if (!empty()) {
31  if(theData.back().recHit()->isValid()) {
32  theNumberOfFoundHits--;
33  theChiSquared -= theData.back().estimate();
34  if(badForCCC(theData.back())) theNumberOfCCCBadHits_--;
35  if(pixel(*(theData.back().recHit()))) theNumberOfFoundPixelHits--;
36  }
37  else if(lost(* (theData.back().recHit()) )) {
38  theNumberOfLostHits--;
39  }
40  else if(isBad(* (theData.back().recHit()) ) && theData.back().recHit()->geographicalId().det()==DetId::Muon ) {
41  theChiSquaredBad -= theData.back().estimate();
42  }
43 
44  theData.pop_back();
45  theNumberOfTrailingFoundHits=countTrailingValidHits(theData);
46  }
47 }
48 
49 
51  push( tm, tm.estimate());
52 }
53 
54 
56  push( tm, tm.estimate());
57 }
58 
59 
60 
61 void Trajectory::push(const TrajectoryMeasurement & tm, double chi2Increment) {
62  theData.push_back(tm); pushAux(chi2Increment);
63 }
64 
65 void Trajectory::push(TrajectoryMeasurement && tm, double chi2Increment) {
66  theData.push_back(tm); pushAux(chi2Increment);
67 }
68 
69 
70 void Trajectory::pushAux(double chi2Increment) {
71  const TrajectoryMeasurement& tm = theData.back();
72  if ( tm.recHit()->isValid()) {
73  theChiSquared += chi2Increment;
74  theNumberOfFoundHits++;
75  theNumberOfTrailingFoundHits++;
76  if (badForCCC(tm)) theNumberOfCCCBadHits_++;
77  if(pixel(*(tm.recHit()))) theNumberOfFoundPixelHits++;
78  }
79  // else if (lost( tm.recHit()) && !inactive(tm.recHit().det())) theNumberOfLostHits++;
80  else if (lost( *(tm.recHit()) ) ) {
81  theNumberOfLostHits++;
82  theNumberOfTrailingFoundHits=0;
83  }
84 
85  else if (isBad( *(tm.recHit()) ) && tm.recHit()->geographicalId().det()==DetId::Muon ) {
86  theChiSquaredBad += chi2Increment;
87  }
88 
89  // in case of a Trajectory constructed without direction,
90  // determine direction from the radii of the first two measurements
91 
92  if ( !theDirectionValidity && theData.size() >= 2) {
93  if (theData[0].updatedState().globalPosition().perp2() <
94  theData.back().updatedState().globalPosition().perp2())
95  theDirection = alongMomentum;
96  else theDirection = oppositeToMomentum;
97  theDirectionValidity = true;
98  }
99 }
100 
101 
102 int Trajectory::ndof(bool bon) const {
103  Trajectory::RecHitContainer && transRecHits = recHits();
104 
105  int dof = 0;
106  int dofBad = 0;
107 
108  for(auto & rechit : transRecHits) {
109  if((rechit)->isValid())
110  dof += (rechit)->dimension();
111  else if( isBad(*rechit) && (rechit)->geographicalId().det()==DetId::Muon )
112  dofBad += (rechit)->dimension();
113  }
114 
115  // If dof!=0 (there is at least 1 valid hit),
116  // return ndof=ndof(fit)
117  // If dof=0 (all rec hits are invalid, only for STA trajectories),
118  // return ndof=ndof(invalid hits)
119  if(dof) {
120  int constr = bon ? 5 : 4;
121  return std::max(dof - constr, 0);
122  }
123  else {
124  // A STA can have < 5 (invalid) hits
125  // if this is the case ==> ndof = 1
126  // (to avoid divisions by 0)
127  int constr = bon ? 5 : 4;
128  return std::max(dofBad - constr, 1);
129  }
130 }
131 
132 
134  hits.reserve(foundHits());
135  for (auto const & tm : theData)
136  if (tm.recHit()->isValid()) hits.push_back(tm.recHit());
137 }
138 
139 
141  if (theDirectionValidity) return theDirection;
142  else throw cms::Exception("TrackingTools/PatternTools","Trajectory::direction() requested but not set");
143 }
144 
145 void Trajectory::check() const {
146  if ( theData.empty())
147  throw cms::Exception("TrackingTools/PatternTools","Trajectory::check() - information requested from empty Trajectory");
148 }
149 
151 {
152  if ( hit.isValid()) return false;
153  else {
154  // // A DetLayer is always inactive in this logic.
155  // // The DetLayer is the Det of an invalid RecHit only if no DetUnit
156  // // is compatible with the predicted state, so we don't really expect
157  // // a hit in this case.
158 
159  if(hit.geographicalId().rawId() == 0) {return false;}
160  else{
161  return hit.getType() == TrackingRecHit::missing;
162  }
163  }
164 }
165 
167 {
168  if ( hit.isValid()) return false;
169  else {
170  if(hit.geographicalId().rawId() == 0) {return false;}
171  else{
172  return hit.getType() == TrackingRecHit::bad;
173  }
174  }
175 }
176 
179  return false;
180  auto const * thit = static_cast<const BaseTrackerRecHit*>( hit.hit() );
181  return thit->isPixel();
182 }
183 
185  if (!trackerHitRTTI::isFromDet(*tm.recHit()))
186  return false;
187  auto const * thit = static_cast<const BaseTrackerRecHit*>( tm.recHit()->hit() );
188  if (!thit)
189  return false;
190  if (thit->isPixel() || thit->isPhase2())
191  return false;
192  if (!tm.updatedState().isValid())
193  return false;
194  return siStripClusterTools::chargePerCM(thit->rawId(),
195  thit->firstClusterRef().stripCluster(),
196  tm.updatedState().localParameters()) < theCCCThreshold_;
197 }
198 
199 void Trajectory::updateBadForCCC(float ccc_threshold) {
200  // If the supplied threshold is the same as the currently cached
201  // one, then return the current number of bad hits for CCC,
202  // otherwise do a new full rescan.
203  if (ccc_threshold == theCCCThreshold_)
204  return;
205 
206  theCCCThreshold_ = ccc_threshold;
207  theNumberOfCCCBadHits_ = 0;
208  for (auto const & h : theData) {
209  if (badForCCC(h))
210  theNumberOfCCCBadHits_++;
211  }
212 }
213 
214 int Trajectory::numberOfCCCBadHits(float ccc_threshold) {
215  updateBadForCCC(ccc_threshold);
216  return theNumberOfCCCBadHits_;
217 }
218 
220 
221  check();
222 
223  //if trajectory is in one half, return the end closer to origin point
224  if ( firstMeasurement().updatedState().globalMomentum().perp() > 1.0
225  && ( firstMeasurement().updatedState().globalPosition().basicVector().dot( firstMeasurement().updatedState().globalMomentum().basicVector() ) *
226  lastMeasurement().updatedState().globalPosition().basicVector().dot( lastMeasurement().updatedState().globalMomentum().basicVector() ) > 0 ) ) {
227  return (firstMeasurement().updatedState().globalPosition().mag() < lastMeasurement().updatedState().globalPosition().mag() ) ?
228  firstMeasurement().updatedState() : lastMeasurement().updatedState();
229  }
230 
231  //more complicated in case of traversing and low-pt trajectories with loops
232  return closestMeasurement(GlobalPoint(0.0,0.0,0.0)).updatedState();
233 
234 }
235 
236 
237 namespace {
239  struct LessMag {
240  LessMag(GlobalPoint point) : thePoint(point) {}
241  bool operator()(const TrajectoryMeasurement& lhs,
242  const TrajectoryMeasurement& rhs) const{
243  if (lhs.updatedState().isValid() && rhs.updatedState().isValid())
244  return (lhs.updatedState().globalPosition() - thePoint).mag2() < (rhs.updatedState().globalPosition() -thePoint).mag2();
245  else
246  {
247  edm::LogError("InvalidStateOnMeasurement")<<"an updated state is not valid. result of LessMag comparator will be wrong.";
248  return false;
249  }
250  }
251  GlobalPoint thePoint;
252  };
253 
254 }
255 
257  check();
258  vector<TrajectoryMeasurement>::const_iterator iter = std::min_element(measurements().begin(), measurements().end(), LessMag(point) );
259 
260  return (*iter);
261 }
262 
264  // reverse the direction (without changing it if it's not along or opposite)
265  if (theDirection == alongMomentum) theDirection = oppositeToMomentum;
266  else if (theDirection == oppositeToMomentum) theDirection = alongMomentum;
267  // reverse the order of the hits
268  std::reverse(theData.begin(), theData.end());
269 }
ConstRecHitPointer const & recHit() const
FWCore Framework interface EventSetupRecordImplementation h
Helper function to determine trigger accepts.
const LocalTrajectoryParameters & localParameters() const
bool isFromDet(TrackingRecHit const &hit)
void pushAux(double chi2Increment)
Definition: Trajectory.cc:70
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
float chargePerCM(DetId detid, Iter a, Iter b)
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
constexpr uint32_t rawId() const
get the raw id
Definition: DetId.h:50
GlobalPoint globalPosition() const
bool badForCCC(const TrajectoryMeasurement &tm)
Definition: Trajectory.cc:184
void check() const
Definition: Trajectory.cc:145
PropagationDirection
TrajectoryMeasurement const & closestMeasurement(GlobalPoint) const
Definition: Trajectory.cc:256
TrajectoryStateOnSurface geometricalInnermostState() const
Definition: Trajectory.cc:219
PropagationDirection const & direction() const
Definition: Trajectory.cc:140
void updateBadForCCC(float ccc_threshold)
Definition: Trajectory.cc:199
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
void validRecHits(ConstRecHitContainer &cont) const
Definition: Trajectory.cc:133
static bool pixel(const TrackingRecHit &hit)
Definition: Trajectory.cc:177
#define end
Definition: vmac.h:39
virtual bool isPixel() const
int ndof(bool bon=true) const
Definition: Trajectory.cc:102
virtual TrackingRecHit const * hit() const
Type getType() const
ConstRecHitContainer RecHitContainer
Definition: Trajectory.h:46
void pop()
Definition: Trajectory.cc:29
static bool isBad(const TrackingRecHit &hit)
Definition: Trajectory.cc:166
bool isValid() const
TrackingRecHit::ConstRecHitContainer ConstRecHitContainer
Definition: Trajectory.h:45
int numberOfCCCBadHits(float ccc_threshold)
Definition: Trajectory.cc:214
T perp() const
Magnitude of transverse component.
#define begin
Definition: vmac.h:32
T dot(const Basic3DVector &v) const
Scalar product, or "dot" product, with a vector of same type.
void reverse()
Definition: Trajectory.cc:263
bool isFromDetOrFast(TrackingRecHit const &hit)
TrajectoryStateOnSurface const & updatedState() const
DetId geographicalId() const
def check(config)
Definition: trackerTree.py:14
uint32_t dimension(pat::CandKinResolution::Parametrization parametrization)
Returns the number of free parameters in a parametrization (3 or 4)
void push(const TrajectoryMeasurement &tm)
Definition: Trajectory.cc:50
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
Definition: invegas.h:5
static bool lost(const TrackingRecHit &hit)
Definition: Trajectory.cc:150