| 1 |
#include "TrackingTools/PatternTools/interface/Trajectory.h"
|
| 2 |
#include "boost/intrusive_ptr.hpp"
|
| 3 |
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
|
| 4 |
#include "FWCore/Utilities/interface/Exception.h"
|
| 5 |
#include "FWCore/MessageLogger/interface/MessageLogger.h"
|
| 6 |
|
| 7 |
#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
|
| 8 |
#include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
|
| 9 |
#include <Geometry/CommonDetUnit/interface/GeomDetUnit.h>
|
| 10 |
#include <Geometry/CommonDetUnit/interface/GeomDetType.h>
|
| 11 |
|
| 12 |
|
| 13 |
#include <algorithm>
|
| 14 |
|
| 15 |
using namespace std;
|
| 16 |
|
| 17 |
void Trajectory::pop() {
|
| 18 |
if (!empty()) {
|
| 19 |
if(theData.back().recHit()->isValid()) {
|
| 20 |
theNumberOfFoundHits--;
|
| 21 |
theChiSquared -= theData.back().estimate();
|
| 22 |
}
|
| 23 |
else if(lost(* (theData.back().recHit()) )) {
|
| 24 |
theNumberOfLostHits--;
|
| 25 |
}
|
| 26 |
else if(isBad(* (theData.back().recHit()) ) && theData.back().recHit()->geographicalId().det()==DetId::Muon ) {
|
| 27 |
theChiSquaredBad -= theData.back().estimate();
|
| 28 |
}
|
| 29 |
|
| 30 |
theData.pop_back();
|
| 31 |
}
|
| 32 |
}
|
| 33 |
|
| 34 |
|
| 35 |
void Trajectory::push( const TrajectoryMeasurement& tm) {
|
| 36 |
push( tm, tm.estimate());
|
| 37 |
}
|
| 38 |
|
| 39 |
void Trajectory::push( const TrajectoryMeasurement& tm, double chi2Increment)
|
| 40 |
{
|
| 41 |
theData.push_back(tm);
|
| 42 |
if ( tm.recHit()->isValid()) {
|
| 43 |
theChiSquared += chi2Increment;
|
| 44 |
theNumberOfFoundHits++;
|
| 45 |
}
|
| 46 |
// else if (lost( tm.recHit()) && !inactive(tm.recHit().det())) theNumberOfLostHits++;
|
| 47 |
else if (lost( *(tm.recHit()) ) ) {
|
| 48 |
theNumberOfLostHits++;
|
| 49 |
}
|
| 50 |
|
| 51 |
else if (isBad( *(tm.recHit()) ) && tm.recHit()->geographicalId().det()==DetId::Muon ) {
|
| 52 |
theChiSquaredBad += chi2Increment;
|
| 53 |
}
|
| 54 |
|
| 55 |
// in case of a Trajectory constructed without direction,
|
| 56 |
// determine direction from the radii of the first two measurements
|
| 57 |
|
| 58 |
if ( !theDirectionValidity && theData.size() >= 2) {
|
| 59 |
if (theData[0].updatedState().globalPosition().perp() <
|
| 60 |
theData.back().updatedState().globalPosition().perp())
|
| 61 |
theDirection = alongMomentum;
|
| 62 |
else theDirection = oppositeToMomentum;
|
| 63 |
theDirectionValidity = true;
|
| 64 |
}
|
| 65 |
}
|
| 66 |
|
| 67 |
Trajectory::RecHitContainer Trajectory::recHits(bool splitting) const {
|
| 68 |
RecHitContainer hits;
|
| 69 |
recHitsV(hits,splitting);
|
| 70 |
return hits;
|
| 71 |
}
|
| 72 |
|
| 73 |
|
| 74 |
int Trajectory::ndof(bool bon) const {
|
| 75 |
const Trajectory::RecHitContainer transRecHits = recHits();
|
| 76 |
|
| 77 |
int dof = 0;
|
| 78 |
int dofBad = 0;
|
| 79 |
|
| 80 |
for(Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin();
|
| 81 |
rechit != transRecHits.end(); ++rechit) {
|
| 82 |
if((*rechit)->isValid())
|
| 83 |
dof += (*rechit)->dimension();
|
| 84 |
else if( isBad(**rechit) && (*rechit)->geographicalId().det()==DetId::Muon )
|
| 85 |
dofBad += (*rechit)->dimension();
|
| 86 |
}
|
| 87 |
|
| 88 |
// If dof!=0 (there is at least 1 valid hit),
|
| 89 |
// return ndof=ndof(fit)
|
| 90 |
// If dof=0 (all rec hits are invalid, only for STA trajectories),
|
| 91 |
// return ndof=ndof(invalid hits)
|
| 92 |
if(dof) {
|
| 93 |
int constr = bon ? 5 : 4;
|
| 94 |
return std::max(dof - constr, 0);
|
| 95 |
}
|
| 96 |
else {
|
| 97 |
// A STA can have < 5 (invalid) hits
|
| 98 |
// if this is the case ==> ndof = 1
|
| 99 |
// (to avoid divisions by 0)
|
| 100 |
int constr = bon ? 5 : 4;
|
| 101 |
return std::max(dofBad - constr, 1);
|
| 102 |
}
|
| 103 |
}
|
| 104 |
|
| 105 |
|
| 106 |
|
| 107 |
void Trajectory::recHitsV(ConstRecHitContainer & hits,bool splitting) const {
|
| 108 |
hits.reserve(theData.size());
|
| 109 |
if(!splitting){
|
| 110 |
for (Trajectory::DataContainer::const_iterator itm
|
| 111 |
= theData.begin(); itm != theData.end(); itm++){
|
| 112 |
hits.push_back((*itm).recHit());
|
| 113 |
}
|
| 114 |
}else{
|
| 115 |
for (Trajectory::DataContainer::const_iterator itm
|
| 116 |
= theData.begin(); itm != theData.end(); itm++){
|
| 117 |
|
| 118 |
// ====== WARNING: this is a temporary solution =========
|
| 119 |
// all this part of code should be implemented internally
|
| 120 |
// in the TrackingRecHit classes. The concrete types of rechit
|
| 121 |
// should be transparent to the Trajectory class
|
| 122 |
|
| 123 |
if( typeid(*(itm->recHit()->hit())) == typeid(SiStripMatchedRecHit2D)){
|
| 124 |
LocalPoint firstLocalPos =
|
| 125 |
itm->updatedState().surface().toLocal(itm->recHit()->transientHits()[0]->globalPosition());
|
| 126 |
|
| 127 |
LocalPoint secondLocalPos =
|
| 128 |
itm->updatedState().surface().toLocal(itm->recHit()->transientHits()[1]->globalPosition());
|
| 129 |
|
| 130 |
LocalVector Delta = secondLocalPos - firstLocalPos;
|
| 131 |
float scalar = Delta.z() * (itm->updatedState().localDirection().z());
|
| 132 |
|
| 133 |
|
| 134 |
TransientTrackingRecHit::ConstRecHitPointer hitA, hitB;
|
| 135 |
|
| 136 |
// Get 2D strip Hits from a matched Hit.
|
| 137 |
//hitA = itm->recHit()->transientHits()[0];
|
| 138 |
//hitB = itm->recHit()->transientHits()[1];
|
| 139 |
|
| 140 |
// Get 2D strip Hits from a matched Hit. Then get the 1D hit from the 2D hit
|
| 141 |
if(!itm->recHit()->transientHits()[0]->detUnit()->type().isEndcap()){
|
| 142 |
hitA = itm->recHit()->transientHits()[0]->transientHits()[0];
|
| 143 |
hitB = itm->recHit()->transientHits()[1]->transientHits()[0];
|
| 144 |
}else{ //don't use 1D hit in the endcap yet
|
| 145 |
hitA = itm->recHit()->transientHits()[0];
|
| 146 |
hitB = itm->recHit()->transientHits()[1];
|
| 147 |
}
|
| 148 |
|
| 149 |
if( (scalar>=0 && direction()==alongMomentum) ||
|
| 150 |
(scalar<0 && direction()==oppositeToMomentum)){
|
| 151 |
hits.push_back(hitA);
|
| 152 |
hits.push_back(hitB);
|
| 153 |
}else if( (scalar>=0 && direction()== oppositeToMomentum) ||
|
| 154 |
(scalar<0 && direction()== alongMomentum)){
|
| 155 |
hits.push_back(hitB);
|
| 156 |
hits.push_back(hitA);
|
| 157 |
}else {
|
| 158 |
//throw cms::Exception("Error in Trajectory::recHitsV(). Direction is not defined");
|
| 159 |
edm::LogError("Trajectory_recHitsV_UndefinedTrackDirection") <<
|
| 160 |
"Error in Trajectory::recHitsV: scalar = " << scalar <<
|
| 161 |
", direction = " << (direction()==alongMomentum ? "along" : (direction()==oppositeToMomentum ? "opposite" : "undefined")) << "\n";
|
| 162 |
hits.push_back(hitA);
|
| 163 |
hits.push_back(hitB);
|
| 164 |
}
|
| 165 |
}else if(typeid(*(itm->recHit()->hit())) == typeid(ProjectedSiStripRecHit2D)){
|
| 166 |
//hits.push_back(itm->recHit()->transientHits()[0]); //Use 2D SiStripRecHit
|
| 167 |
if(!itm->recHit()->transientHits()[0]->detUnit()->type().isEndcap()){
|
| 168 |
hits.push_back(itm->recHit()->transientHits()[0]->transientHits()[0]); //Use 1D SiStripRecHit
|
| 169 |
}else{
|
| 170 |
hits.push_back(itm->recHit()->transientHits()[0]); //Use 2D SiStripRecHit
|
| 171 |
}
|
| 172 |
// ===================================================================================
|
| 173 |
}else if(typeid(*(itm->recHit()->hit())) == typeid(SiStripRecHit2D)){
|
| 174 |
//hits.push_back(itm->recHit()); //Use 2D SiStripRecHit
|
| 175 |
if(!itm->recHit()->detUnit()->type().isEndcap()){
|
| 176 |
hits.push_back(itm->recHit()->transientHits()[0]); //Use 1D SiStripRecHit
|
| 177 |
}else{
|
| 178 |
hits.push_back(itm->recHit()); //Use 2D SiStripRecHit
|
| 179 |
}
|
| 180 |
}else{
|
| 181 |
hits.push_back(itm->recHit());
|
| 182 |
}
|
| 183 |
}//end loop on measurements
|
| 184 |
}
|
| 185 |
}
|
| 186 |
|
| 187 |
void Trajectory::validRecHits(ConstRecHitContainer & hits) const {
|
| 188 |
hits.reserve(foundHits());
|
| 189 |
for (Trajectory::DataContainer::const_iterator itm
|
| 190 |
= theData.begin(); itm != theData.end(); itm++)
|
| 191 |
if ((*itm).recHit()->isValid()) hits.push_back((*itm).recHit());
|
| 192 |
}
|
| 193 |
|
| 194 |
|
| 195 |
PropagationDirection const & Trajectory::direction() const {
|
| 196 |
if (theDirectionValidity) return theDirection;
|
| 197 |
else throw cms::Exception("TrackingTools/PatternTools","Trajectory::direction() requested but not set");
|
| 198 |
}
|
| 199 |
|
| 200 |
void Trajectory::check() const {
|
| 201 |
if ( theData.empty())
|
| 202 |
throw cms::Exception("TrackingTools/PatternTools","Trajectory::check() - information requested from empty Trajectory");
|
| 203 |
}
|
| 204 |
|
| 205 |
bool Trajectory::lost( const TransientTrackingRecHit& hit)
|
| 206 |
{
|
| 207 |
if ( hit.isValid()) return false;
|
| 208 |
else {
|
| 209 |
// // A DetLayer is always inactive in this logic.
|
| 210 |
// // The DetLayer is the Det of an invalid RecHit only if no DetUnit
|
| 211 |
// // is compatible with the predicted state, so we don't really expect
|
| 212 |
// // a hit in this case.
|
| 213 |
|
| 214 |
if(hit.geographicalId().rawId() == 0) {return false;}
|
| 215 |
else{
|
| 216 |
return hit.getType() == TrackingRecHit::missing;
|
| 217 |
}
|
| 218 |
}
|
| 219 |
}
|
| 220 |
|
| 221 |
bool Trajectory::isBad( const TransientTrackingRecHit& hit)
|
| 222 |
{
|
| 223 |
if ( hit.isValid()) return false;
|
| 224 |
else {
|
| 225 |
if(hit.geographicalId().rawId() == 0) {return false;}
|
| 226 |
else{
|
| 227 |
return hit.getType() == TrackingRecHit::bad;
|
| 228 |
}
|
| 229 |
}
|
| 230 |
}
|
| 231 |
|
| 232 |
TrajectoryStateOnSurface Trajectory::geometricalInnermostState() const {
|
| 233 |
|
| 234 |
check();
|
| 235 |
|
| 236 |
//if trajectory is in one half, return the end closer to origin point
|
| 237 |
if ( firstMeasurement().updatedState().globalMomentum().perp() > 1.0
|
| 238 |
&& ( firstMeasurement().updatedState().globalPosition().basicVector().dot( firstMeasurement().updatedState().globalMomentum().basicVector() ) *
|
| 239 |
lastMeasurement().updatedState().globalPosition().basicVector().dot( lastMeasurement().updatedState().globalMomentum().basicVector() ) > 0 ) ) {
|
| 240 |
return (firstMeasurement().updatedState().globalPosition().mag() < lastMeasurement().updatedState().globalPosition().mag() ) ?
|
| 241 |
firstMeasurement().updatedState() : lastMeasurement().updatedState();
|
| 242 |
}
|
| 243 |
|
| 244 |
//more complicated in case of traversing and low-pt trajectories with loops
|
| 245 |
return closestMeasurement(GlobalPoint(0.0,0.0,0.0)).updatedState();
|
| 246 |
|
| 247 |
}
|
| 248 |
|
| 249 |
|
| 250 |
namespace {
|
| 251 |
/// used to determine closest measurement to given point
|
| 252 |
struct LessMag {
|
| 253 |
LessMag(GlobalPoint point) : thePoint(point) {}
|
| 254 |
bool operator()(const TrajectoryMeasurement& lhs,
|
| 255 |
const TrajectoryMeasurement& rhs) const{
|
| 256 |
if (lhs.updatedState().isValid() && rhs.updatedState().isValid())
|
| 257 |
return (lhs.updatedState().globalPosition() - thePoint).mag2() < (rhs.updatedState().globalPosition() -thePoint).mag2();
|
| 258 |
else
|
| 259 |
{
|
| 260 |
edm::LogError("InvalidStateOnMeasurement")<<"an updated state is not valid. result of LessMag comparator will be wrong.";
|
| 261 |
return false;
|
| 262 |
}
|
| 263 |
}
|
| 264 |
GlobalPoint thePoint;
|
| 265 |
};
|
| 266 |
|
| 267 |
}
|
| 268 |
|
| 269 |
TrajectoryMeasurement const & Trajectory::closestMeasurement(GlobalPoint point) const {
|
| 270 |
check();
|
| 271 |
vector<TrajectoryMeasurement>::const_iterator iter = std::min_element(measurements().begin(), measurements().end(), LessMag(point) );
|
| 272 |
|
| 273 |
return (*iter);
|
| 274 |
}
|
| 275 |
|
| 276 |
void Trajectory::reverse() {
|
| 277 |
// reverse the direction (without changing it if it's not along or opposite)
|
| 278 |
if (theDirection == alongMomentum) theDirection = oppositeToMomentum;
|
| 279 |
else if (theDirection == oppositeToMomentum) theDirection = alongMomentum;
|
| 280 |
// reverse the order of the hits
|
| 281 |
std::reverse(theData.begin(), theData.end());
|
| 282 |
}
|
