4 added + 1 modified, total 5 files
trf++/src/gtrfit
diff -N BiGTrackFitter.cpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ BiGTrackFitter.cpp 15 Aug 2011 17:09:48 -0000 1.1
@@ -0,0 +1,231 @@
+// BiGTrackFitter.cpp
+
+#include "gtrfit/BiGTrackFitter.hpp"
+#include <cassert>
+#include <cmath>
+#include <map>
+#include "trfutil/STLprint.h"
+#include "trfbase/SurfacePtr.h"
+#include "trfdca/SurfDCA.h"
+
+using std::cout;
+using std::endl;
+using std::map;
+using std::abs;
+using trf::PropagatorPtr;
+using trf::SurfacePtr;
+using trf::SurfDCA;
+using trf::MissPtr;
+
+//**********************************************************************
+
+// Constructor.
+
+BiGTrackFitter::BiGTrackFitter(const PropagatorPtr& pprop,
+//cng D0ClusterRegistry* preg,
+ double chsq_max,
+ int maxloop)
+: _infit(pprop, -1, /*cng preg, */ chsq_max),
+ _outfit(pprop, 1, /*cng preg,*/ chsq_max),
+ _maxloop(maxloop)
+ { }
+
+
+//**********************************************************************
+
+// Backward compatibility constructor.
+
+BiGTrackFitter::BiGTrackFitter(const PropagatorPtr& pprop
+//cng ,D0ClusterRegistry* preg)
+ )
+: _infit(pprop, -1, /*cng preg,*/ 100.),
+ _outfit(pprop, 1, /*cng preg,*/ 100.),
+ _maxloop(1)
+ { }
+
+
+//**********************************************************************
+
+// Fitter.
+
+int BiGTrackFitter::fit(GTrack& trg) const {
+
+ // Loop over the GTrack states and count how many states have valid
+ // clusters
+ int nstates=0;
+ const GTrack::StateList& states = trg.get_states();
+ GTrack::StateList::const_iterator istate;
+ for( istate = states.begin() ; istate != states.end() ; ++istate )
+ if( istate->cluster() )
+ nstates++;
+ if(nstates < 2)
+ return 500;
+
+ int nloop = _maxloop;
+ if ( nloop <= 0 ) {
+ nloop = 1;
+ }
+ GTrack trg_in = trg;
+ GTrack trg_out;
+ for ( int iloop=0; iloop<nloop; ++iloop ) {
+
+ // Fit inside to out.
+ trg_out = trg_in;
+ {
+ int stat = _outfit.fit(trg_out);
+ if ( stat ) return 200 + stat;
+ if ( trg_out.get_states().size() != trg.get_states().size() ) {
+ return 302;
+ }
+ }
+
+ // Fit outside to in.
+ trg_in = trg_out;
+ {
+ int stat = _infit.fit(trg_in);
+ if ( stat ) return 100 + stat;
+ if ( trg_in.get_states().size() != trg.get_states().size() ) {
+ return 301;
+ }
+ }
+
+
+ }
+
+ // Fetch the state lists.
+ const GTrack::StateList& out_states = trg_out.get_states();
+ map<double, GTrackState> in_state_map;
+ {
+ GTrack::StateList in_states = trg_in.get_states();
+ assert( in_states.size() == out_states.size() );
+
+ // Match states from the in and out fits. We do this by matching the
+ // surface, the cluster index, and the miss, if any. Move the in_states
+ // into a map indexed by the path distance from the out_states.
+
+ for(GTrack::StateList::const_iterator i = out_states.begin();
+ i != out_states.end(); ++i) {
+ assert(i->has_valid_fit());
+ double s = i->s();
+//cng const ChunkClusterIndex& ccidx = i->cluster();
+ const trf::ClusterPtr pclu = i->cluster();
+ const MissPtr& pmiss = i->miss();
+ const SurfacePtr& psrf = i->track().get_surface();
+ bool at_dca = psrf->get_type() == SurfDCA::get_static_type();
+//cng if(ccidx.is_valid()) {
+ if(pclu) {
+ assert(pmiss == 0);
+ assert(psrf->get_type() != SurfDCA::get_static_type());
+ }
+ else {
+ if(pmiss == 0)
+ assert(at_dca);
+ else
+ assert(!at_dca);
+ }
+ bool match = false;
+ for(GTrack::StateList::iterator j = in_states.begin();
+ j != in_states.end(); ++j) {
+ assert(j->has_valid_fit());
+ match =
+//cng ccidx.is_valid() && ccidx == j->cluster() ||
+ pclu && pclu == j->cluster() ||
+ pmiss != 0 && j->miss() != 0 && *pmiss == *(j->miss()) ||
+ at_dca && j->track().get_surface()->get_type() ==
+ SurfDCA::get_static_type();
+ if(match) {
+ in_state_map[s] = *j;
+ in_states.erase(j);
+ break;
+ }
+ }
+ assert(match);
+ }
+ assert(in_state_map.size() == out_states.size());
+ assert(in_states.size() == 0);
+ }
+
+ // Build the new state list.
+ GTrack::StateList newstates;
+
+ // Create states using s from the inside->out track and the fit
+ // from the outside-to-in track until we find a cluster. We
+ // expect the fits will all be optimal.
+ //
+ GTrack::StateList::const_iterator iout_state = out_states.begin();
+ while ( iout_state != out_states.end() ) {
+ const GTrackState& outstate = *iout_state;
+ const GTrackState& instate = in_state_map[outstate.s()];
+ ++iout_state;
+ GTrackState state;
+//cng if(instate.cluster().is_valid()) {
+ if(instate.cluster()) {
+ assert(instate.miss() == 0);
+ state = GTrackState(outstate.s(), instate.track(), instate.fit_status(),
+ instate.chi_square(), instate.cluster());
+ }
+ else
+ state = GTrackState(outstate.s(), instate.track(), instate.fit_status(),
+ instate.chi_square(), instate.miss());
+ if ( state.fit_status() != GTrackState::OPTIMAL ) {
+ cout << endl << "BiGTrackFitter: Nonoptimal state." << endl;
+ cout << endl;
+ print_list("Out states", out_states);
+ cout << endl;
+ cout << endl << "BiGTrackFitter: End error message." << endl;
+ assert(false);
+ }
+ assert( state.fit_status() == GTrackState::OPTIMAL );
+ newstates.insert(state);
+ // Check if this is the first outside-to-in cluster.
+ // If so, check the same is true fo inside-to-out and then exit
+ // the loop.
+ if ( instate.cluster() ) {
+ // Make sure inside-to-out tracks has the same cluster.
+ if ( ! outstate.cluster() ) {
+ return 303;
+ }
+ assert( instate.cluster() == outstate.cluster() );
+ break;
+ // If outside-to-in is not a cluster, verify the same is true
+ // for inside-to-out.
+ } else {
+ assert( outstate.miss() != 0 ||
+ outstate.track().get_surface()->get_type() ==
+ SurfDCA::get_static_type());
+ if ( instate.cluster() ) {
+ return 304;
+ }
+ assert( instate.miss() != 0 ||
+ instate.track().get_surface()->get_type() ==
+ SurfDCA::get_static_type());
+ }
+ }
+ if ( iout_state == out_states.end() ) {
+ cout << endl << "BiGTrackFitter: Bad state selection." << endl;
+ cout << endl << out_states << endl;
+ cout << endl << "BiGTrackFitter: End error message." << endl;
+ }
+ assert( iout_state != out_states.end() );
+
+ // Take the remaining states from the inside->out fit.
+ while ( iout_state != out_states.end() ) {
+ const GTrackState& outstate = *iout_state;
+ const GTrackState& instate = in_state_map[outstate.s()];
+ ++iout_state;
+ newstates.insert(outstate);
+ }
+ assert( iout_state == out_states.end() );
+
+ assert( newstates.size() == trg.get_states().size() );
+ assert( newstates.begin()->fit_status() == GTrackState::OPTIMAL );
+ assert( newstates.rbegin()->fit_status() == GTrackState::OPTIMAL );
+
+ double bfield = trg_in.get_bfield();
+ trg = GTrack(newstates, bfield);
+//cng trg.set_quality(trg_in.get_quality());
+ return 0;
+
+}
+
+//**********************************************************************
trf++/src/gtrfit
diff -N GTrackResidual.cpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ GTrackResidual.cpp 15 Aug 2011 17:09:48 -0000 1.1
@@ -0,0 +1,315 @@
+// GTrackResidual.cpp
+
+#include "gtrfit/GTrackResidual.hpp"
+#include <cassert>
+#include "trfbase/Surface.h"
+#include "trfbase/Hit.h"
+#include "trfbase/ETrack.h"
+#include "trfbase/PropStat.h"
+#include "trfbase/Propagator.h"
+#include "trfbase/PropagatorPtr.h"
+#include "trfbase/Hit.h"
+#include "gtrbase/GTrackState.hpp"
+#include "gtrbase/GTrack.hpp"
+#include "gtrfit/SimpleGTrackFitter.hpp"
+//cng #include "d0cluster_evt/ChunkClusterPtr.hpp"
+
+using trf::TrackVector;
+using trf::TrackError;
+using trf::TrackDerivative;
+using trf::Surface;
+using trf::SurfacePtr;
+using trf::ETrack;
+using trf::Hit;
+using trf::HitList;
+using trf::ClusterPtr;
+using trf::PropStat;
+using trf::Propagator;
+using trf::PropagatorPtr;
+
+//**********************************************************************
+// Local definitions.
+//**********************************************************************
+
+namespace {
+
+//**********************************************************************
+
+// Maximum chi-square for refit.
+double csqmax = 1000000.0;
+
+//**********************************************************************
+
+// Implementation class.
+
+class GTrackResidualImp {
+
+public: // data
+
+ PropagatorPtr _pprop;
+ ETrack* _ptre;
+ double _chisq;
+ HitList _hits;
+ SimpleGTrackFitter _ffit;
+ SimpleGTrackFitter _bfit;
+
+public: // functions
+
+ GTrackResidualImp(const PropagatorPtr& pprop) :
+ _pprop(pprop), _ptre(0), _chisq(-1.),
+ _ffit(_pprop, 1, csqmax),
+ _bfit(_pprop, -1, csqmax) { }
+
+
+ ~GTrackResidualImp() {
+ delete _ptre;
+ }
+
+};
+
+//**********************************************************************
+
+} // end unnamed namespace.
+
+//**********************************************************************
+// Member functions.
+//**********************************************************************
+
+GTrackResidual::GTrackResidual(const PropagatorPtr& pprop)
+: _pimp(new GTrackResidualImp(pprop)) {
+ _pimp->_ptre = 0;
+}
+
+//**********************************************************************
+
+// Destructor.
+
+GTrackResidual::~GTrackResidual() {
+ delete _pimp;
+}
+
+//**********************************************************************
+
+// Calculate a residual.
+
+int GTrackResidual::calc(const GTrack& trg, const GTrackState& rstate) {
+
+ ETrack*& ptre = _pimp->_ptre;
+ double& chisq = _pimp->_chisq;
+ HitList& hits = _pimp->_hits;
+ const Propagator& prop = *_pimp->_pprop;
+ const SimpleGTrackFitter& ffit = _pimp->_ffit;
+ const SimpleGTrackFitter& bfit = _pimp->_bfit;
+
+ // Clean up from last call.
+ if ( ptre != 0 ) {
+ delete ptre;
+ ptre = 0;
+ }
+ chisq = -1.;
+
+ // Fetch the list of states.
+ GTrack::StateList states = trg.get_states();
+
+ // Build list of states in front of the residual state.
+ GTrack::StateList fstates;
+ GTrack::StateList::const_iterator istate = states.begin();
+ while ( istate != states.end() ) {
+ if ( *istate == rstate ) {
+ break;
+ }
+ fstates.insert(*istate);
+ ++istate;
+ }
+
+ // Exit if state was not found.
+ if ( istate == states.end () ) {
+ return 1;
+ }
+
+ // Build the list of states after the residual state.
+ GTrack::StateList bstates;
+ while ( ++istate != states.end() ) {
+ if ( *istate == rstate ) {
+ return 2;
+ }
+ bstates.insert(*istate);
+ }
+
+ assert( fstates.size()+bstates.size()+1 == states.size() );
+
+ // Build and refit the front track.
+ ETrack ftre;
+ double fchisq = -1.;
+ int nfclu = 0;
+ for( GTrack::StateList::const_iterator f = fstates.begin();
+ f != fstates.end(); ++f)
+ if(f->cluster())
+ ++nfclu;
+ if ( nfclu ) {
+ GTrack ftrg(fstates);
+ int stat = ffit.fit(ftrg);
+ if ( stat != 0 ) {
+ return 1000 + stat;
+ }
+ GTrackState fstate = *ftrg.get_states().rbegin();
+ assert( fstate.fit_status() == GTrackState::OPTIMAL );
+ ftre = fstate.track();
+ fchisq = fstate.chi_square();
+ }
+
+ // Build and refit the back track.
+ ETrack btre;
+ double bchisq = -1.;
+ int nbclu = 0;
+ for( GTrack::StateList::const_iterator b = bstates.begin();
+ b != bstates.end(); ++b)
+ if(b->cluster())
+ ++nbclu;
+ if ( nbclu ) {
+ GTrack btrg(bstates);
+ int stat = bfit.fit(btrg);
+ if ( stat != 0 ) {
+ return 2000 + stat;
+ }
+ GTrackState bstate = *btrg.get_states().begin();
+ assert( bstate.fit_status() == GTrackState::OPTIMAL );
+ btre = bstate.track();
+ bchisq = bstate.chi_square();
+ }
+
+ // Extract the TRF++ cluster from the residual state.
+//cng const ChunkClusterPtr pd0clu = rstate.cluster();
+//cng ClusterPtr pclu(pd0clu->new_trf_cluster());
+
+ const ClusterPtr pclu = rstate.cluster();
+ // Propagate the tracks to the cluster surface.
+ const Surface& srf = pclu->get_surface();
+ if ( ftre.is_valid() ) {
+ PropStat pstat = prop.err_dir_prop(ftre, srf, Propagator::FORWARD);
+ if ( ! pstat.success() ) {
+ return 3;
+ }
+ }
+ if ( btre.is_valid() ) {
+ PropStat pstat = prop.err_dir_prop(btre, srf, Propagator::BACKWARD);
+ if ( ! pstat.success() ) {
+ return 4;
+ }
+ }
+
+
+ // Average the tracks.
+ // First state.
+ if ( ! btre.is_valid() ) {
+ assert( ftre.is_valid() );
+ ptre = new ETrack(ftre);
+ chisq = fchisq;
+
+ // Last state.
+ } else if ( ! ftre.is_valid() ) {
+ assert( btre.is_valid() );
+ ptre = new ETrack(btre);
+ chisq = bchisq;
+
+ // Intermediate state.
+ } else {
+ // Fetch the track vectors.
+ const TrackVector& fvec = ftre.get_vector();
+ const TrackVector& bvec = btre.get_vector();
+ // Fetch the error matrices.
+ const TrackError& ferr = ftre.get_error();
+ const TrackError& berr = btre.get_error();
+
+ // Difference vector.
+ TrackVector dvec = fvec - bvec;
+
+ // Inverse of error sum
+ TrackError tinv = ferr + berr;
+ {
+ int istat = geninvert(tinv);
+ if ( istat != 0 ) {
+ return 3000 + istat;
+ }
+ }
+
+ // Calculate new chisquare.
+ chisq = fchisq + bchisq + tinv % dvec;
+
+ // Define average track vector and error matrix;
+ TrackVector avec;
+ TrackError aerr;
+
+ // Calculate perturbation on better-measured track.
+ TrackVector tdvec = tinv * dvec;
+ if(fstates.size() > bstates.size()) {
+ TrackVector ftdvec = ferr * tdvec;
+ avec = fvec - ftdvec;
+ TrackDerivative ferr0(ferr);
+ TrackError ftf(ferr0 % tinv);
+ aerr = ferr - ftf;
+ }
+ else { // fstates.size() <= bstates.size()
+ TrackVector btdvec = berr * tdvec;
+ avec = bvec + btdvec;
+ TrackDerivative berr0(berr);
+ TrackError btb(berr0 % tinv);
+ aerr = berr - btb;
+ }
+
+ // Build the new track.
+ SurfacePtr psrf(srf.new_pure_surface());
+ ptre = new ETrack(psrf, avec, aerr);
+ // Some surfaces require that we set the direction.
+ if ( ! ptre->is_valid() ) {
+ if ( ftre.is_forward() ) {
+ assert( btre.is_forward() );
+ ptre->set_forward();
+ } else if ( ftre.is_backward() ) {
+ assert( btre.is_backward() );
+ ptre->set_backward();
+ } else {
+ assert(false);
+ }
+ }
+ assert( ptre->is_valid() );
+ }
+ assert( ptre!=0 && ptre->is_valid() );
+
+ // Predict a hit with the average track.
+ hits = pclu->predict(*ptre);
+ if ( hits.size() != 1 ) {
+ return 6;
+ }
+
+ return 0;
+
+}
+
+//**********************************************************************
+
+// Return the track.
+
+const ETrack& GTrackResidual::track() const {
+ assert( _pimp->_ptre != 0 );
+ return *_pimp->_ptre;
+}
+
+//**********************************************************************
+
+// Return the chisquare.
+
+double GTrackResidual::chisq() const {
+ return _pimp->_chisq;
+}
+
+//**********************************************************************
+
+// Return the hit.
+
+const Hit& GTrackResidual::hit() const {
+ assert( _pimp->_hits.size() == 1 );
+ return *_pimp->_hits.front();
+}
+
+//**********************************************************************
trf++/src/gtrfit
diff -N GTrackSmoothingRefitter.cpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ GTrackSmoothingRefitter.cpp 15 Aug 2011 17:09:48 -0000 1.1
@@ -0,0 +1,254 @@
+#include "gtrfit/GTrackSmoothingRefitter.hpp"
+//#include "gtrbase/GtrTrackQuality.hpp"
+#include "trfdca/SurfDCA.h"
+#include <cstdlib>
+
+using namespace trf;
+
+//********************************************************************
+
+namespace {
+ const GTrackState& find_state(const GTrackState& state,const GTrack& trgs0,const GTrack& trgs1);
+ const GTrackState& find_state(const GTrackState& state,const GTrack& trgs0);
+ GTrack reorder_gtrack(const GTrack& gtr0,const GTrack& trgs0);
+ GTrack reorder(const GTrack& gtr0);
+}
+
+//**************************************************************************
+
+GTrackSmoothingRefitter::GTrackSmoothingRefitter(const trf::PropagatorPtr& pprop, int direction, double chsq ) : _fitter_inout(pprop,1,chsq),_fitter_outin(pprop,-1,chsq),_smoother(pprop) {
+}
+
+//**************************************************************************
+
+int GTrackSmoothingRefitter::fit(GTrack& trg) const {
+
+ GTrack trg0 = trg;
+
+ // Loop over the GTrack states and count how many states have valid
+ // clusters
+ int nstates=0;
+ const GTrack::StateList& states = trg.get_states();
+ GTrack::StateList::const_iterator istate;
+ for( istate = states.begin() ; istate != states.end() ; ++istate )
+ if( istate->cluster() )
+ nstates++;
+
+ int half0 = nstates-2;
+ if(nstates >= 6)
+ half0 = nstates-3;
+ if(nstates >= 8)
+ half0 = nstates-4;
+ int half1 = nstates-half0;
+ if(half0 <= 0 || half1 <= 0)
+ return 500;
+
+ int stat = _fitter_inout.fit(trg0);
+
+ if( stat )
+ return 200 + stat;
+
+ trg0 = reorder(trg0);
+
+ GTrack trg0s = trg0;
+
+ stat = _smoother.smooth(trg0s,half0);
+
+ if( stat )
+ return 400 + stat;
+
+ GTrack trg10s = trg0;
+
+ stat = _fitter_outin.fit(trg10s);
+
+ if( stat )
+ return 100 + stat;
+
+ // reorder trg
+ GTrack trg1s = reorder_gtrack(trg0,trg10s);
+ stat = _smoother.smooth(trg1s,half1);
+
+ if( stat )
+ return 300 + stat;
+
+ // combined trg0s and trg1s (different smoothed parts of the same track
+ // into a single track using s from trg0 (INOUT fit)
+ const GTrack::StateList& states0 = trg0.get_states();
+
+ GTrack::StateList newstates;
+
+ for( istate = states0.begin() ; istate != states0.end() ; ++istate ) {
+ const GTrackState& state = find_state(*istate,trg0s,trg1s);
+
+ GTrackState newstate;
+
+ if( state.cluster()) {
+ assert( state.miss() == 0);
+ newstate = GTrackState(state, istate->s());
+ }
+ else
+ newstate = GTrackState(state, istate->s());
+
+ newstates.insert(newstate);
+
+ }
+
+
+ // assert( newstates.size() == trg0s.get_states().size() + trg1s.get_states().size() );
+
+ // Make a copy of the quality from the original track.
+ //cng GtrTrackQuality quality = trg.get_quality();
+
+ // Construct new GTrack.
+ double bfield = trg.get_bfield();
+ trg = GTrack(newstates,bfield);
+
+ // Reinsert quality object.
+ //cng trg.set_quality(quality);
+
+ return stat;
+
+}
+
+//********************************************************************
+
+namespace {
+ const GTrackState& find_state(const GTrackState& state,const GTrack& trgs0,const GTrack& trgs1) {
+
+ const GTrack::StateList& out_states = trgs0.get_states();
+ const GTrack::StateList& in_states = trgs1.get_states();
+
+
+ assert(state.has_valid_fit());
+ double s = state.s();
+// const ChunkClusterIndex& ccidx = state.cluster();
+ const trf::ClusterPtr& pclu = state.cluster();
+ const MissPtr& pmiss = state.miss();
+ bool at_dca =
+ state.track().get_surface()->get_type() == SurfDCA::get_static_type();
+ if(pclu)
+ assert(pmiss == 0 && !at_dca);
+ else {
+ if(pmiss == 0)
+ assert(at_dca);
+ else
+ assert(!at_dca);
+ }
+ bool match = false;
+
+ for(GTrack::StateList::const_iterator j = in_states.begin();
+ j != in_states.end() ; ++j) {
+ assert(j->has_valid_fit());
+ match =
+ pclu && pclu == j->cluster() ||
+ pmiss != 0 && j->miss() != 0 && *pmiss == *(j->miss()) ||
+ at_dca && j->track().get_surface()->get_type() ==
+ SurfDCA::get_static_type();
+ if( match )
+ return *j;
+ }
+
+ for(GTrack::StateList::const_iterator j = out_states.begin();
+ j != out_states.end() ; ++j) {
+ assert(j->has_valid_fit());
+ match =
+ pclu && pclu == j->cluster() ||
+ pmiss != 0 && j->miss() != 0 && *pmiss == *(j->miss()) ||
+ at_dca && j->track().get_surface()->get_type() ==
+ SurfDCA::get_static_type();
+ if( match )
+ return *j;
+ }
+
+ assert(false);
+ static GTrackState invalid;
+ return invalid;
+ }
+
+const GTrackState& find_state(const GTrackState& state,const GTrack& trgs0) {
+
+ const GTrack::StateList& in_states = trgs0.get_states();
+
+ assert(state.has_valid_fit());
+ double s = state.s();
+
+//cng const ChunkClusterIndex& ccidx = state.cluster();
+ const trf::ClusterPtr pclu = state.cluster();
+ const MissPtr& pmiss = state.miss();
+ bool at_dca =
+ state.track().get_surface()->get_type() == SurfDCA::get_static_type();
+ if(pclu)
+ assert(pmiss == 0 && !at_dca);
+ else {
+ if(pmiss == 0)
+ assert(at_dca);
+ else
+ assert(!at_dca);
+ }
+ bool match = false;
+
+ for(GTrack::StateList::const_iterator j = in_states.begin();
+ j != in_states.end() ; ++j) {
+ assert(j->has_valid_fit());
+ match =
+//cng ccidx.is_valid() && ccidx == j->cluster() ||
+ pclu && pclu == j->cluster() ||
+ pmiss != 0 && j->miss() != 0 && *pmiss == *(j->miss()) ||
+ at_dca && j->track().get_surface()->get_type() ==
+ SurfDCA::get_static_type();
+ if( match )
+ return *j;
+ }
+
+ assert(false);
+ static GTrackState invalid;
+ return invalid;
+ }
+
+GTrack reorder_gtrack(const GTrack& gtr0,const GTrack& trgs0) {
+ const GTrack::StateList& states = gtr0.get_states();
+ GTrack::StateList newstates;
+ for(GTrack::StateList::const_iterator i=states.begin() ; i!=states.end(); i++ ) {
+ const GTrackState& state = find_state(*i,trgs0);
+
+ GTrackState newstate;
+ if( state.cluster()) {
+ assert( state.miss() == 0);
+ newstate = GTrackState(state, i->s());
+ }
+ else
+ newstate = GTrackState(state, i->s());
+
+ newstates.insert(newstate);
+ }
+ double bfield = gtr0.get_bfield();
+ return GTrack(newstates,bfield);
+}
+
+GTrack reorder(const GTrack& gtr0) {
+ const GTrack::StateList& states = gtr0.get_states();
+ GTrack::StateList newstates;
+ double s_old=0.;
+ for(GTrack::StateList::const_iterator i=states.begin() ; i!=states.end(); i++ ) {
+ const GTrackState& state = *i;
+ double s = state.s();
+ if( abs(s-s_old) < 1e-4 && i!=states.begin() )
+ s+=2e-4;
+ s_old = s;
+ GTrackState newstate;
+ if( state.cluster() ) {
+ assert( state.miss() == 0);
+ newstate = GTrackState(state, s);
+ }
+ else
+ newstate = GTrackState(state, s);
+
+ newstates.insert(newstate);
+ }
+
+ double bfield = gtr0.get_bfield();
+ return GTrack(newstates,bfield);
+}
+
+}
+//********************************************************************
trf++/src/gtrfit
diff -N GTrackSmoothingRefitter2.cpp
--- /dev/null 1 Jan 1970 00:00:00 -0000
+++ GTrackSmoothingRefitter2.cpp 15 Aug 2011 17:09:48 -0000 1.1
@@ -0,0 +1,340 @@
+#include "gtrfit/GTrackSmoothingRefitter2.hpp"
+#include "trfdca/SurfDCA.h"
+
+using namespace trf;
+
+//********************************************************************
+
+namespace {
+ const GTrackState& find_state(const GTrackState& state, const GTrack& trgs0);
+ GTrackState merge_states(const GTrackState&, const GTrackState&);
+ GTrack reorder_gtrack(const GTrack& gtr0, const GTrack& trgs0);
+ GTrack reorder(const GTrack& gtr0);
+ GTrack merge_gtracks(const GTrack&, const GTrack&);
+}
+
+//**************************************************************************
+
+GTrackSmoothingRefitter2::
+GTrackSmoothingRefitter2(const trf::PropagatorPtr& pprop,
+ /*cng D0ClusterRegistry* reg,*/ double chsq) :
+ _fitter_inout(pprop,1,/*cng reg,*/chsq),
+ _fitter_outin(pprop,-1,/*cng reg,*/ chsq)
+{}
+
+//**************************************************************************
+
+int GTrackSmoothingRefitter2::fit(GTrack& trg) const {
+
+ // Do in->out fit.
+
+ GTrack trg0 = trg;
+ int stat = _fitter_inout.fit(trg0);
+
+ if( stat )
+ return 200 + stat;
+
+ // Make sure each state has a unique distance s with a minimum separation.
+
+ trg0 = reorder(trg0);
+
+ // Do out->in fit.
+
+ GTrack trg1 = trg0;
+ stat = _fitter_outin.fit(trg1);
+
+ if( stat )
+ return 100 + stat;
+
+ // Reorder trg1 so that states are in the same order as trg0.
+
+ trg1 = reorder_gtrack(trg0, trg1);
+
+ // Combine states from trg0 and trg1.
+
+ trg = merge_gtracks(trg0, trg1);
+
+ if(!trg.is_valid())
+ return 300;
+
+ return stat;
+}
+
+//********************************************************************
+
+namespace {
+
+ //
+ // This function locates the state from a GTrack which matches a given
+ // state (same surface, same cluster index, same miss).
+ //
+
+ const GTrackState& find_state(const GTrackState& state, const GTrack& trg) {
+
+ const GTrack::StateList& in_states = trg.get_states();
+
+ assert(state.has_valid_fit());
+
+ // Extract information from input state.
+
+ double s = state.s();
+//cng const ChunkClusterIndex& ccidx = state.cluster();
+ const trf::ClusterPtr& pclu = state.cluster();
+ const MissPtr& pmiss = state.miss();
+ bool at_dca =
+ state.track().get_surface()->get_type() == SurfDCA::get_static_type();
+
+ // Only one of cluster, miss, dca.
+
+ if(pclu)
+ assert(pmiss == 0 && !at_dca);
+ else {
+ if(pmiss == 0)
+ assert(at_dca);
+ else
+ assert(!at_dca);
+ }
+
+ // Search for matching state in GTrack.
+
+ bool match = false;
+ for(GTrack::StateList::const_iterator j = in_states.begin();
+ j != in_states.end() ; ++j) {
+ assert(j->has_valid_fit());
+ match =
+//cng ccidx.is_valid() && ccidx == j->cluster() ||
+ pclu && pclu == j->cluster() ||
+ pmiss != 0 && j->miss() != 0 && *pmiss == *(j->miss()) ||
+ at_dca && j->track().get_surface()->get_type() ==
+ SurfDCA::get_static_type();
+ if( match )
+ return *j;
+ }
+
+ // Should never exit the loop.
+
+ assert(false);
+ static GTrackState invalid;
+ return invalid;
+ }
+
+
+ //
+ // This function creates a new GTrack, taking the distances from the first
+ // argument, and tracks and other state attributes from corresponding
+ // surfaces from the second argument.
+ //
+
+ GTrack reorder_gtrack(const GTrack& trg0, const GTrack& trg1) {
+
+ GTrack::StateList newstates;
+
+ // Loop over states from first GTrack. These states supply the distance
+ // and define the order.
+
+ const GTrack::StateList& states0 = trg0.get_states();
+ for(GTrack::StateList::const_iterator i=states0.begin(); i!=states0.end();
+ ++i ) {
+
+ // Find corresponding state from second track.
+
+ const GTrackState& state = find_state(*i, trg1);
+
+ // Make a new state, taking distance from first track, and other
+ // attributes from second track.
+
+ newstates.insert(GTrackState(state, i->s()));
+ }
+ double bfield = trg0.get_bfield();
+ GTrack new_trg(newstates,bfield);
+//cng new_trg.set_quality(trg0.get_quality());
+ return new_trg;
+ }
+
+
+ //
+ // This function merges (averages) two states. In case of error, a
+ // default-constructed state is returned. A default-constructed state
+ // has fit_status BADSTATE, which can be tested using method is_valid().
+ //
+
+ GTrackState merge_states(const GTrackState& state0,
+ const GTrackState& state1) {
+
+ // If one of the two states is already optimal, just return that state.
+
+ GTrackState::FitStatus stat0 = state0.fit_status();
+ GTrackState::FitStatus stat1 = state1.fit_status();
+ if(stat0 == GTrackState::OPTIMAL)
+ return state0;
+ if(stat1 == GTrackState::OPTIMAL)
+ return state1;
+
+ // Otherwise, the only legitimate states are forward + backward.
+
+ assert(stat0 == GTrackState::FORWARD && stat1 == GTrackState::BACKWARD ||
+ stat0 == GTrackState::BACKWARD && stat1 == GTrackState::FORWARD);
+
+ // Extract linear algebra information.
+
+ const ETrack& tre0 = state0.track();
+ const ETrack& tre1 = state1.track();
+ const TrackVector *pveca, *pvecb;
+ const TrackError *perra, *perrb;
+ double chia, chib;
+
+ // Track "a" is the better measured track.
+ // We use cached no-hit information from the less well measured track.
+
+ if(tre0.get_error(4,4) < tre1.get_error(4,4)) {
+ pveca = &tre0.get_vector();
+ perra = &tre0.get_error();
+ chia = state0.chi_square();
+ assert(!state1.cluster() || state1.has_nohit_track());
+ const ETrack& tre1_nohit = state1.get_nohit_track();
+ pvecb = &tre1_nohit.get_vector();
+ perrb = &tre1_nohit.get_error();
+ chib = state1.get_nohit_chi_square();
+ }
+ else {
+ pveca = &tre1.get_vector();
+ perra = &tre1.get_error();
+ chia = state1.chi_square();
+ assert(!state0.cluster() || state0.has_nohit_track());
+ const ETrack& tre0_nohit = state0.get_nohit_track();
+ pvecb = &tre0_nohit.get_vector();
+ perrb = &tre0_nohit.get_error();
+ chib = state0.get_nohit_chi_square();
+ }
+
+ // Calculate average track quantities.
+
+ const Surface& srf0 = *tre0.get_surface();
+ TrackVector diff = srf0.vec_diff(*pveca, *pvecb); // veca-vecb w/phi wrap
+ TrackError tinv = *perra + *perrb;
+ int istat = geninvert(tinv);
+ if(istat)
+ return GTrackState();
+
+ // Average chisquare.
+
+ double chi = chia + chib + tinv % diff;
+
+ // Average track parameters.
+
+ TrackVector dvec = *perra * (tinv * diff);
+ TrackVector vec = *pveca - dvec;
+
+ // Average track error.
+
+ TrackDerivative erra(*perra); // Non-symmetric representation.
+ TrackError derr(erra % tinv);
+ TrackError err = *perra - derr;
+
+ // Make new ETrack.
+
+ SurfacePtr psrf(srf0.new_pure_surface());
+ ETrack tre(psrf, vec, err);
+ if (!tre.is_valid()) {
+ if(tre0.is_forward()) {
+ if(!tre1.is_forward())
+ return GTrackState();
+ tre.set_forward();
+ }
+ else if(tre0.is_backward()) {
+ if(!tre1.is_backward())
+ return GTrackState();
+ tre.set_backward();
+ }
+ }
+ assert(tre.is_valid());
+
+ // Now put everything together into a new state.
+
+//cng const ChunkClusterIndex& icclu = state0.cluster();
+ const trf::ClusterPtr& pclu = state0.cluster();
+ MissPtr pmiss = state0.miss();
+//cng bool has_cluster = icclu.is_valid();
+ bool has_cluster = pclu;
+ bool has_miss = pmiss != 0;
+
+ // Cluster state.
+
+ if(has_cluster) {
+ assert(pclu == state1.cluster());
+ return GTrackState(state0.s(), tre, GTrackState::OPTIMAL, chi, pclu);
+ }
+
+ // Miss state.
+
+ else if(has_miss) {
+ assert(*pmiss == *state1.miss());
+ return GTrackState(state0.s(), tre, GTrackState::OPTIMAL, chi, pmiss);
+ }
+
+ // Default (no cluster, no miss -- can occur with dca state).
+
+ return GTrackState(state0.s(), tre, GTrackState::OPTIMAL, chi);
+ }
+
+
+ //
+ // This function checks that adjacent states have a certain minimum
+ // separation in distance s.
+ //
+
+ GTrack reorder(const GTrack& trg) {
+ const GTrack::StateList& states = trg.get_states();
+ GTrack::StateList newstates;
+ double s_old=0.;
+ for(GTrack::StateList::const_iterator i=states.begin() ; i!=states.end();
+ ++i ) {
+ const GTrackState& state = *i;
+ double s = state.s();
+ if( i!=states.begin()) {
+ double smin = s_old + 1.e-4;
+ if(s < smin)
+ s = smin;
+ }
+ s_old = s;
+ newstates.insert(GTrackState(state, s));
+ }
+
+ double bfield = trg.get_bfield();
+ GTrack new_trg(newstates,bfield);
+//cng new_trg.set_quality(trg.get_quality());
+ return new_trg;
+ }
+
+
+ //
+ // This function combines two GTracks that were fit in opposite directions.
+ // The states must have already been reordered so that corresponding
+ // states are in the same order between the two tracks.
+ //
+
+ GTrack merge_gtracks(const GTrack& trg0, const GTrack& trg1) {
+
+ GTrack::StateList newstates;
+
+ const GTrack::StateList states0 = trg0.get_states();
+ const GTrack::StateList states1 = trg1.get_states();
+ GTrack::StateList::const_iterator i = states0.begin();
+ GTrack::StateList::const_iterator j = states1.begin();
+
+ for(; i != states0.end(); ++i, ++j ) {
+ assert(i->s() == j->s());
+ GTrackState state = merge_states(*i, *j);
+ if(!state.is_valid())
+ return GTrack();
+ newstates.insert(state);
+ }
+
+ double bfield = trg0.get_bfield();
+ GTrack new_trg(newstates,bfield);
+//cng new_trg.set_quality(trg0.get_quality());
+ return new_trg;
+ }
+}
+
+//********************************************************************
trf++/src/gtrfit
diff -u -r1.1 -r1.2
--- COMPONENTS 11 Aug 2011 00:04:40 -0000 1.1
+++ COMPONENTS 15 Aug 2011 17:09:48 -0000 1.2
@@ -1,2 +1,6 @@
GTrackSmoother
SimpleGTrackFitter
+GTrackSmoothingRefitter
+GTrackSmoothingRefitter2
+BiGTrackFitter
+GTrackResidual
CVSspam 0.2.8