PatternTools/src/Trajectory.cc

#include "TrackingTools/PatternTools/interface/Trajectory.h"
#include "boost/intrusive_ptr.hpp" 
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h" 
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2D.h"
#include "DataFormats/TrackerRecHit2D/interface/ProjectedSiStripRecHit2D.h"
#include <Geometry/CommonDetUnit/interface/GeomDetUnit.h>
#include <Geometry/CommonDetUnit/interface/GeomDetType.h>


#include <algorithm>

using namespace std;

void Trajectory::pop() {
  if (!empty()) {
    if(theData.back().recHit()->isValid()) {
      theNumberOfFoundHits--;
      theChiSquared -= theData.back().estimate();
    }
    else if(lost(* (theData.back().recHit()) )) {
      theNumberOfLostHits--;
    }
    else if(isBad(* (theData.back().recHit()) ) && theData.back().recHit()->geographicalId().det()==DetId::Muon ) {
      theChiSquaredBad -= theData.back().estimate();
    }

    theData.pop_back();
  }
}


void Trajectory::push( const TrajectoryMeasurement& tm) {
  push( tm, tm.estimate());
}

void Trajectory::push( const TrajectoryMeasurement& tm, double chi2Increment)
{
  theData.push_back(tm);
  if ( tm.recHit()->isValid()) {
    theChiSquared += chi2Increment;
    theNumberOfFoundHits++;
  }
  // else if (lost( tm.recHit()) && !inactive(tm.recHit().det())) theNumberOfLostHits++;
  else if (lost( *(tm.recHit()) ) ) {
    theNumberOfLostHits++;
  }
 
  else if (isBad( *(tm.recHit()) ) && tm.recHit()->geographicalId().det()==DetId::Muon ) {
    theChiSquaredBad += chi2Increment;
  }
 
  // in case of a Trajectory constructed without direction, 
  // determine direction from the radii of the first two measurements

  if ( !theDirectionValidity && theData.size() >= 2) {
    if (theData[0].updatedState().globalPosition().perp() <
        theData.back().updatedState().globalPosition().perp())
      theDirection = alongMomentum;
    else theDirection = oppositeToMomentum;
    theDirectionValidity = true;
  }
}

Trajectory::RecHitContainer Trajectory::recHits(bool splitting) const {
  RecHitContainer hits;
  recHitsV(hits,splitting);
  return hits;
}


int Trajectory::ndof(bool bon) const {
  const Trajectory::RecHitContainer transRecHits = recHits();
  
  int dof = 0;
  int dofBad = 0;
  
  for(Trajectory::RecHitContainer::const_iterator rechit = transRecHits.begin();
      rechit != transRecHits.end(); ++rechit) {
    if((*rechit)->isValid())
      dof += (*rechit)->dimension();
    else if( isBad(**rechit) && (*rechit)->geographicalId().det()==DetId::Muon )
      dofBad += (*rechit)->dimension();
  }

  // If dof!=0 (there is at least 1 valid hit),
  //    return ndof=ndof(fit)
  // If dof=0 (all rec hits are invalid, only for STA trajectories),
  //    return ndof=ndof(invalid hits)
  if(dof) {
    int constr = bon ? 5 : 4;
    return std::max(dof - constr, 0);
  }
  else {
    // A STA can have < 5 (invalid) hits
    // if this is the case ==> ndof = 1
    // (to avoid divisions by 0)
    int constr = bon ? 5 : 4;
    return std::max(dofBad - constr, 1);
  }
}


void Trajectory::recHitsV(ConstRecHitContainer & hits,bool splitting) const {
  hits.reserve(theData.size());
  if(!splitting){  
    for (Trajectory::DataContainer::const_iterator itm
           = theData.begin(); itm != theData.end(); itm++){    
      hits.push_back((*itm).recHit());
    }
  }else{    
    for (Trajectory::DataContainer::const_iterator itm
           = theData.begin(); itm != theData.end(); itm++){    

      // ====== WARNING: this is a temporary solution =========
      //        all this part of code should be implemented internally 
      //        in the TrackingRecHit classes. The concrete types of rechit 
      //        should be transparent to the Trajectory class

      if( typeid(*(itm->recHit()->hit())) == typeid(SiStripMatchedRecHit2D)){
        LocalPoint firstLocalPos = 
          itm->updatedState().surface().toLocal(itm->recHit()->transientHits()[0]->globalPosition());
        
        LocalPoint secondLocalPos = 
          itm->updatedState().surface().toLocal(itm->recHit()->transientHits()[1]->globalPosition());
        
        LocalVector Delta = secondLocalPos - firstLocalPos;
        float scalar  = Delta.z() * (itm->updatedState().localDirection().z());
        

        TransientTrackingRecHit::ConstRecHitPointer hitA, hitB;

        // Get 2D strip Hits from a matched Hit.
        //hitA = itm->recHit()->transientHits()[0];
        //hitB = itm->recHit()->transientHits()[1];

        // Get 2D strip Hits from a matched Hit. Then get the 1D hit from the 2D hit
        if(!itm->recHit()->transientHits()[0]->detUnit()->type().isEndcap()){
          hitA = itm->recHit()->transientHits()[0]->transientHits()[0];
          hitB = itm->recHit()->transientHits()[1]->transientHits()[0];
        }else{ //don't use 1D hit in the endcap yet
          hitA = itm->recHit()->transientHits()[0];
          hitB = itm->recHit()->transientHits()[1];
        }

        if( (scalar>=0 && direction()==alongMomentum) ||
            (scalar<0 && direction()==oppositeToMomentum)){
          hits.push_back(hitA);
          hits.push_back(hitB);
        }else if( (scalar>=0 && direction()== oppositeToMomentum) ||
                  (scalar<0 && direction()== alongMomentum)){
          hits.push_back(hitB);
          hits.push_back(hitA);
        }else {
          //throw cms::Exception("Error in Trajectory::recHitsV(). Direction is not defined");  
          edm::LogError("Trajectory_recHitsV_UndefinedTrackDirection") << 
                "Error in Trajectory::recHitsV: scalar = " << scalar << 
                ", direction = " << (direction()==alongMomentum ? "along" : (direction()==oppositeToMomentum ? "opposite" : "undefined")) << "\n";
          hits.push_back(hitA);
          hits.push_back(hitB);
        }         
      }else if(typeid(*(itm->recHit()->hit())) == typeid(ProjectedSiStripRecHit2D)){
        //hits.push_back(itm->recHit()->transientHits()[0]);    //Use 2D SiStripRecHit
        if(!itm->recHit()->transientHits()[0]->detUnit()->type().isEndcap()){
          hits.push_back(itm->recHit()->transientHits()[0]->transientHits()[0]);        //Use 1D SiStripRecHit
        }else{
          hits.push_back(itm->recHit()->transientHits()[0]);    //Use 2D SiStripRecHit
        }
        // ===================================================================================  
      }else if(typeid(*(itm->recHit()->hit())) == typeid(SiStripRecHit2D)){
        //hits.push_back(itm->recHit());  //Use 2D SiStripRecHit
        if(!itm->recHit()->detUnit()->type().isEndcap()){
          hits.push_back(itm->recHit()->transientHits()[0]); //Use 1D SiStripRecHit
        }else{
          hits.push_back(itm->recHit());  //Use 2D SiStripRecHit
        }
      }else{
        hits.push_back(itm->recHit());
      }
    }//end loop on measurements
  }
}

void Trajectory::validRecHits(ConstRecHitContainer & hits) const {
  hits.reserve(foundHits());
  for (Trajectory::DataContainer::const_iterator itm
         = theData.begin(); itm != theData.end(); itm++)
    if ((*itm).recHit()->isValid()) hits.push_back((*itm).recHit());
}


PropagationDirection const & Trajectory::direction() const {
  if (theDirectionValidity) return theDirection;
  else throw cms::Exception("TrackingTools/PatternTools","Trajectory::direction() requested but not set");
}

void Trajectory::check() const {
  if ( theData.empty()) 
    throw cms::Exception("TrackingTools/PatternTools","Trajectory::check() - information requested from empty Trajectory");
}

bool Trajectory::lost( const TransientTrackingRecHit& hit)
{
  if ( hit.isValid()) return false;
  else {
  //     // A DetLayer is always inactive in this logic.
  //     // The DetLayer is the Det of an invalid RecHit only if no DetUnit 
  //     // is compatible with the predicted state, so we don't really expect
  //     // a hit in this case.
  
    if(hit.geographicalId().rawId() == 0) {return false;}
    else{
      return hit.getType() == TrackingRecHit::missing;
    }
  }
}

bool Trajectory::isBad( const TransientTrackingRecHit& hit)
{
  if ( hit.isValid()) return false;
  else {
    if(hit.geographicalId().rawId() == 0) {return false;}
    else{
      return hit.getType() == TrackingRecHit::bad;
    }
  }
}

TrajectoryStateOnSurface Trajectory::geometricalInnermostState() const {

  check();

  //if trajectory is in one half, return the end closer to origin point
  if ( firstMeasurement().updatedState().globalMomentum().perp() > 1.0
      && ( firstMeasurement().updatedState().globalPosition().basicVector().dot( firstMeasurement().updatedState().globalMomentum().basicVector() ) *
       lastMeasurement().updatedState().globalPosition().basicVector().dot( lastMeasurement().updatedState().globalMomentum().basicVector() )  > 0 ) ) {
     return (firstMeasurement().updatedState().globalPosition().mag() < lastMeasurement().updatedState().globalPosition().mag() ) ?
            firstMeasurement().updatedState() : lastMeasurement().updatedState();
  }

  //more complicated in case of traversing and low-pt trajectories with loops
  return closestMeasurement(GlobalPoint(0.0,0.0,0.0)).updatedState();

}


namespace {
  /// used to determine closest measurement to given point
  struct LessMag {
    LessMag(GlobalPoint point) : thePoint(point) {}
    bool operator()(const TrajectoryMeasurement& lhs,
                    const TrajectoryMeasurement& rhs) const{ 
      if (lhs.updatedState().isValid() && rhs.updatedState().isValid())
        return (lhs.updatedState().globalPosition() - thePoint).mag2() < (rhs.updatedState().globalPosition() -thePoint).mag2();
      else
        {
          edm::LogError("InvalidStateOnMeasurement")<<"an updated state is not valid. result of LessMag comparator will be wrong.";
          return false;
        }
    }
    GlobalPoint thePoint;
  };

}

TrajectoryMeasurement const & Trajectory::closestMeasurement(GlobalPoint point) const {
  check();
  vector<TrajectoryMeasurement>::const_iterator iter = std::min_element(measurements().begin(), measurements().end(), LessMag(point) );

  return (*iter);
}

void Trajectory::reverse() {
    // reverse the direction (without changing it if it's not along or opposite)
    if (theDirection == alongMomentum)           theDirection = oppositeToMomentum;
    else if (theDirection == oppositeToMomentum) theDirection = alongMomentum;
    // reverse the order of the hits
    std::reverse(theData.begin(), theData.end());
}
Revision:	1.1
Committed:	Fri Nov 25 16:38:25 2011 UTC (13 years, 5 months ago) by econte
Content type:	text/plain
Branch:	MAIN
CVS Tags:	TBD2011, TBD_2011, HEAD
Log Message:	new IPHC alignment
#	User	Rev	Content
1	econte	1.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			}