Track.cc

/* Copyright 2008-2010, Technische Universitaet Muenchen,
   Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
 
   This file is part of GENFIT.
 
   GENFIT is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published
   by the Free Software Foundation, either version 3 of the License, or
   (at your option) any later version.
 
   GENFIT is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.
 
   You should have received a copy of the GNU Lesser General Public License
   along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
*/
 
#include "Track.h"
 
#include "Exception.h"
#include "IO.h"
#include "KalmanFitterInfo.h"
#include "KalmanFitStatus.h"
#include "PlanarMeasurement.h"
#include "AbsMeasurement.h"
 
#include "WireTrackCandHit.h"
 
#include <algorithm>
#include <map>
 
#include <TDatabasePDG.h>
#include <TMath.h>
#include <TBuffer.h>
 
//#include <glog/logging.h>
 
//#define DEBUG
 
 
namespace genfit {
 
  Track::Track() :
    cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(6), covSeed_(6)
  {
    ;
  }
 
 
  Track::Track(const TrackCand& trackCand, const MeasurementFactory<AbsMeasurement>& factory, AbsTrackRep* rep) :
    cardinalRep_(0), fitStatuses_(), stateSeed_(6), covSeed_(6)
  {
 
    if (rep != nullptr)
      addTrackRep(rep);
 
    createMeasurements(trackCand, factory);
 
    // Copy seed information from candidate
    timeSeed_ = trackCand.getTimeSeed();
    stateSeed_ = trackCand.getStateSeed();
    covSeed_ = trackCand.getCovSeed();
 
    mcTrackId_ = trackCand.getMcTrackId();
 
    // fill cache
    fillPointsWithMeasurement();
 
    checkConsistency();
  }
 
  void
  Track::createMeasurements(const TrackCand& trackCand, const MeasurementFactory<AbsMeasurement>& factory)
  {
    // create the measurements using the factory.
    std::vector <AbsMeasurement*> factoryHits = factory.createMany(trackCand);
 
    if (factoryHits.size() != trackCand.getNHits()) {
      Exception exc("Track::Track ==> factoryHits.size() != trackCand->getNHits()", __LINE__, __FILE__);
      exc.setFatal();
      throw exc;
    }
 
    // create TrackPoints
    for (unsigned int i = 0; i < factoryHits.size(); ++i) {
      TrackPoint* tp = new TrackPoint(factoryHits[i], this);
      tp->setSortingParameter(trackCand.getHit(i)->getSortingParameter());
      insertPoint(tp);
    }
  }
 
 
  Track::Track(AbsTrackRep* trackRep, const TVectorD& stateSeed) :
    cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(stateSeed),
    covSeed_(TMatrixDSym::kUnit, TMatrixDSym(6))
  {
    addTrackRep(trackRep);
  }
 
 
  Track::Track(AbsTrackRep* trackRep, const TVector3& posSeed, const TVector3& momSeed) :
    cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(6),
    covSeed_(TMatrixDSym::kUnit, TMatrixDSym(6))
  {
    addTrackRep(trackRep);
    setStateSeed(posSeed, momSeed);
  }
 
 
  Track::Track(AbsTrackRep* trackRep, const TVectorD& stateSeed, const TMatrixDSym& covSeed) :
    cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(stateSeed),
    covSeed_(covSeed)
  {
    addTrackRep(trackRep);
  }
 
 
  Track::Track(const Track& rhs) :
    TObject(rhs),
    cardinalRep_(rhs.cardinalRep_), mcTrackId_(rhs.mcTrackId_), timeSeed_(rhs.timeSeed_),
    stateSeed_(rhs.stateSeed_), covSeed_(rhs.covSeed_)
  {
    rhs.checkConsistency();
 
    std::map<const AbsTrackRep*, AbsTrackRep*> oldRepNewRep;
 
    for (std::vector<AbsTrackRep*>::const_iterator it = rhs.trackReps_.begin(); it != rhs.trackReps_.end(); ++it) {
      AbsTrackRep* newRep = (*it)->clone();
      addTrackRep(newRep);
      oldRepNewRep[(*it)] = newRep;
    }
 
    trackPoints_.reserve(rhs.trackPoints_.size());
    for (std::vector<TrackPoint*>::const_iterator it = rhs.trackPoints_.begin(); it != rhs.trackPoints_.end(); ++it) {
      trackPoints_.push_back(new TrackPoint(**it, oldRepNewRep));
      trackPoints_.back()->setTrack(this);
    }
 
    for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it = rhs.fitStatuses_.begin(); it != rhs.fitStatuses_.end(); ++it) {
      setFitStatus(it->second->clone(), oldRepNewRep[it->first]);
    }
 
    fillPointsWithMeasurement();
 
    checkConsistency();
  }
 
  Track& Track::operator=(Track other)
  {
    swap(other);
 
    for (std::vector<TrackPoint*>::const_iterator it = trackPoints_.begin(); it != trackPoints_.end(); ++it) {
      (*it)->setTrack(this);
    }
 
    fillPointsWithMeasurement();
 
    checkConsistency();
 
    return *this;
  }
 
  void Track::swap(Track& other)
  {
    std::swap(this->trackReps_, other.trackReps_);
    std::swap(this->cardinalRep_, other.cardinalRep_);
    std::swap(this->trackPoints_, other.trackPoints_);
    std::swap(this->trackPointsWithMeasurement_, other.trackPointsWithMeasurement_);
    std::swap(this->fitStatuses_, other.fitStatuses_);
    std::swap(this->mcTrackId_, other.mcTrackId_);
    std::swap(this->timeSeed_, other.timeSeed_);
    std::swap(this->stateSeed_, other.stateSeed_);
    std::swap(this->covSeed_, other.covSeed_);
 
  }
 
  Track::~Track()
  {
    this->Clear();
  }
 
  void Track::Clear(Option_t*)
  {
    // This function is needed for TClonesArray embedding.
    // FIXME: smarter containers or pointers needed ...
    for (size_t i = 0; i < trackPoints_.size(); ++i)
      delete trackPoints_[i];
 
    trackPoints_.clear();
    trackPointsWithMeasurement_.clear();
 
    for (std::map< const AbsTrackRep*, FitStatus* >::iterator it = fitStatuses_.begin(); it != fitStatuses_.end(); ++it)
      delete it->second;
    fitStatuses_.clear();
 
    for (size_t i = 0; i < trackReps_.size(); ++i)
      delete trackReps_[i];
    trackReps_.clear();
 
    cardinalRep_ = 0;
 
    mcTrackId_ = -1;
 
    timeSeed_ = 0;
    stateSeed_.Zero();
    covSeed_.Zero();
  }
 
 
  TrackPoint* Track::getPoint(int id) const
  {
    if (id < 0)
      id += trackPoints_.size();
 
    return trackPoints_.at(id);
  }
 
 
  TrackPoint* Track::getPointWithMeasurement(int id) const
  {
    if (id < 0)
      id += trackPointsWithMeasurement_.size();
 
    return trackPointsWithMeasurement_.at(id);
  }
 
 
  TrackPoint* Track::getPointWithMeasurementAndFitterInfo(int id, const AbsTrackRep* rep) const
  {
    if (rep == nullptr)
      rep = getCardinalRep();
 
    if (id >= 0) {
      int i = 0;
      for (std::vector<TrackPoint*>::const_iterator it = trackPointsWithMeasurement_.begin(); it != trackPointsWithMeasurement_.end();
           ++it) {
        if ((*it)->hasFitterInfo(rep)) {
          if (id == i)
            return (*it);
          ++i;
        }
      }
    } else {
      // Search backwards.
      int i = -1;
      for (std::vector<TrackPoint*>::const_reverse_iterator it = trackPointsWithMeasurement_.rbegin();
           it != trackPointsWithMeasurement_.rend(); ++it) {
        if ((*it)->hasFitterInfo(rep)) {
          if (id == i)
            return (*it);
          --i;
        }
      }
    }
 
    // Not found, i.e. abs(id) > number of fitted TrackPoints
    return 0;
  }
 
 
  TrackPoint* Track::getPointWithFitterInfo(int id, const AbsTrackRep* rep) const
  {
    if (rep == nullptr)
      rep = getCardinalRep();
 
    if (id >= 0) {
      int i = 0;
      for (std::vector<TrackPoint*>::const_iterator it = trackPoints_.begin(); it != trackPoints_.end(); ++it) {
        if ((*it)->hasFitterInfo(rep)) {
          if (id == i)
            return (*it);
          ++i;
        }
      }
    } else {
      // Search backwards.
      int i = -1;
      for (std::vector<TrackPoint*>::const_reverse_iterator it = trackPoints_.rbegin(); it != trackPoints_.rend(); ++it) {
        if ((*it)->hasFitterInfo(rep)) {
          if (id == i)
            return (*it);
          --i;
        }
      }
    }
 
    // Not found, i.e. abs(id) > number of fitted TrackPoints
    return 0;
  }
 
 
  const MeasuredStateOnPlane& Track::getFittedState(int id, const AbsTrackRep* rep, bool biased) const
  {
    if (rep == nullptr)
      rep = getCardinalRep();
 
    TrackPoint* point = getPointWithFitterInfo(id, rep);
    if (point == nullptr) {
      Exception exc("Track::getFittedState ==> no trackPoint with fitterInfo for rep", __LINE__, __FILE__);
      exc.setFatal();
      throw exc;
    }
    return point->getFitterInfo(rep)->getFittedState(biased);
  }
 
 
  int Track::getIdForRep(const AbsTrackRep* rep) const
  {
    for (size_t i = 0; i < trackReps_.size(); ++i)
      if (trackReps_[i] == rep)
        return i;
 
    Exception exc("Track::getIdForRep ==> cannot find TrackRep in Track", __LINE__, __FILE__);
    exc.setFatal();
    throw exc;
  }
 
 
  bool Track::hasFitStatus(const AbsTrackRep* rep) const
  {
    if (rep == nullptr)
      rep = getCardinalRep();
 
    if (fitStatuses_.find(rep) == fitStatuses_.end())
      return false;
 
    return (fitStatuses_.at(rep) != nullptr);
  }
 
 
  bool Track::hasKalmanFitStatus(const AbsTrackRep* rep) const
  {
    if (rep == nullptr)
      rep = getCardinalRep();
 
    if (fitStatuses_.find(rep) == fitStatuses_.end())
      return false;
 
    return (dynamic_cast<KalmanFitStatus*>(fitStatuses_.at(rep)) != nullptr);
  }
 
 
  KalmanFitStatus* Track::getKalmanFitStatus(const AbsTrackRep* rep) const
  {
    return dynamic_cast<KalmanFitStatus*>(getFitStatus(rep));
  }
 
 
  void Track::setFitStatus(FitStatus* fitStatus, const AbsTrackRep* rep)
  {
    if (fitStatuses_.find(rep) != fitStatuses_.end())
      delete fitStatuses_.at(rep);
 
    fitStatuses_[rep] = fitStatus;
  }
 
 
  void Track::setStateSeed(const TVector3& pos, const TVector3& mom)
  {
    stateSeed_.ResizeTo(6);
 
    stateSeed_(0) = pos.X();
    stateSeed_(1) = pos.Y();
    stateSeed_(2) = pos.Z();
 
    stateSeed_(3) = mom.X();
    stateSeed_(4) = mom.Y();
    stateSeed_(5) = mom.Z();
  }
 
 
 
  void Track::insertPoint(TrackPoint* point, int id)
  {
 
    point->setTrack(this);
 
#ifdef DEBUG
    debugOut << "Track::insertPoint at position " << id  << "\n";
#endif
    assert(point != nullptr);
    trackHasChanged();
 
    point->setTrack(this);
 
    if (trackPoints_.size() == 0) {
      trackPoints_.push_back(point);
 
      if (point->hasRawMeasurements())
        trackPointsWithMeasurement_.push_back(point);
 
      return;
    }
 
    if (id == -1 || id == (int)trackPoints_.size()) {
      trackPoints_.push_back(point);
 
      if (point->hasRawMeasurements())
        trackPointsWithMeasurement_.push_back(point);
 
      deleteReferenceInfo(std::max(0, (int)trackPoints_.size() - 2), (int)trackPoints_.size() - 1);
 
      // delete fitter infos if inserted point has a measurement
      if (point->hasRawMeasurements()) {
        deleteForwardInfo(-1, -1);
        deleteBackwardInfo(0, -2);
      }
 
      return;
    }
 
    // [-size, size-1] is the allowed range
    assert(id < (ssize_t)trackPoints_.size() || -id - 1 <= (ssize_t)trackPoints_.size());
 
    if (id < 0)
      id += trackPoints_.size() + 1;
 
    // insert
    trackPoints_.insert(trackPoints_.begin() + id, point);  // insert inserts BEFORE
 
    // delete fitter infos if inserted point has a measurement
    if (point->hasRawMeasurements()) {
      deleteForwardInfo(id, -1);
      deleteBackwardInfo(0, id);
    }
 
    // delete reference info of neighbouring points
    deleteReferenceInfo(std::max(0, id - 1), std::min((int)trackPoints_.size() - 1, id + 1));
 
    fillPointsWithMeasurement();
  }
 
 
  void Track::insertPoints(std::vector<TrackPoint*> points, int id)
  {
 
    int nBefore = getNumPoints();
    int n = points.size();
 
    if (n == 0)
      return;
    if (n == 1) {
      insertPoint(points[0], id);
      return;
    }
 
    for (std::vector<TrackPoint*>::iterator p = points.begin(); p != points.end(); ++p)
      (*p)->setTrack(this);
 
    if (id == -1 || id == (int)trackPoints_.size()) {
      trackPoints_.insert(trackPoints_.end(), points.begin(), points.end());
 
      deleteReferenceInfo(std::max(0, nBefore - 1), nBefore);
 
      deleteForwardInfo(nBefore, -1);
      deleteBackwardInfo(0, std::max(0, nBefore - 1));
 
      fillPointsWithMeasurement();
 
      return;
    }
 
 
    assert(id < (ssize_t)trackPoints_.size() || -id - 1 <= (ssize_t)trackPoints_.size());
 
    if (id < 0)
      id += trackPoints_.size() + 1;
 
 
    // insert
    trackPoints_.insert(trackPoints_.begin() + id, points.begin(), points.end());  // insert inserts BEFORE
 
    // delete fitter infos if inserted point has a measurement
    deleteForwardInfo(id, -1);
    deleteBackwardInfo(0, id + n);
 
    // delete reference info of neighbouring points
    deleteReferenceInfo(std::max(0, id - 1), std::min((int)trackPoints_.size() - 1, id));
    deleteReferenceInfo(std::max(0, id + n - 1), std::min((int)trackPoints_.size() - 1, id + n));
 
    fillPointsWithMeasurement();
  }
 
 
  void Track::deletePoint(int id)
  {
 
#ifdef DEBUG
    debugOut << "Track::deletePoint at position " << id  << "\n";
#endif
 
    trackHasChanged();
 
    if (id < 0)
      id += trackPoints_.size();
    assert(id > 0);
 
 
    // delete forwardInfo after point (backwardInfo before point) if deleted point has a measurement
    if (trackPoints_[id]->hasRawMeasurements()) {
      deleteForwardInfo(id, -1);
      deleteBackwardInfo(0, id - 1);
    }
 
    // delete reference info of neighbouring points
    deleteReferenceInfo(std::max(0, id - 1), std::min((int)trackPoints_.size() - 1, id + 1));
 
    // delete point
    std::vector<TrackPoint*>::iterator it = std::find(trackPointsWithMeasurement_.begin(), trackPointsWithMeasurement_.end(),
                                                      trackPoints_[id]);
    if (it != trackPointsWithMeasurement_.end())
      trackPointsWithMeasurement_.erase(it);
 
    delete trackPoints_[id];
    trackPoints_.erase(trackPoints_.begin() + id);
 
    fillPointsWithMeasurement();
 
  }
 
 
  void Track::insertMeasurement(AbsMeasurement* measurement, int id)
  {
    insertPoint(new TrackPoint(measurement, this), id);
  }
 
  void Track::deleteFittedState(const genfit::AbsTrackRep* rep)
  {
    if (hasFitStatus(rep)) {
      delete fitStatuses_.at(rep);
      fitStatuses_.erase(rep);
    }
 
    // delete FitterInfos related to the deleted TrackRep
    for (const auto& trackPoint : trackPoints_) {
      if (trackPoint->hasFitterInfo(rep)) {
        trackPoint->deleteFitterInfo(rep);
      }
    }
  }
 
 
  void Track::mergeTrack(const Track* other, int id)
  {
 
#ifdef DEBUG
    debugOut << "Track::mergeTrack\n";
#endif
 
    if (other->getNumPoints() == 0)
      return;
 
    std::map<const AbsTrackRep*, AbsTrackRep*> otherRepThisRep;
    std::vector<const AbsTrackRep*> otherRepsToRemove;
 
    for (std::vector<AbsTrackRep*>::const_iterator otherRep = other->trackReps_.begin(); otherRep != other->trackReps_.end();
         ++otherRep) {
      bool found(false);
      for (std::vector<AbsTrackRep*>::const_iterator thisRep = trackReps_.begin(); thisRep != trackReps_.end(); ++thisRep) {
        if ((*thisRep)->isSame(*otherRep)) {
          otherRepThisRep[*otherRep] = *thisRep;
#ifdef DEBUG
          debugOut << " map other rep " << *otherRep << " to " << (*thisRep) << "\n";
#endif
          if (found) {
            Exception exc("Track::mergeTrack ==> more than one matching rep.", __LINE__, __FILE__);
            exc.setFatal();
            throw exc;
          }
          found = true;
          break;
        }
      }
      if (!found) {
        otherRepsToRemove.push_back(*otherRep);
#ifdef DEBUG
        debugOut << " remove other rep " << *otherRep << "\n";
#endif
      }
    }
 
 
    std::vector<TrackPoint*> points;
    points.reserve(other->getNumPoints());
 
    for (std::vector<TrackPoint*>::const_iterator otherTp = other->trackPoints_.begin(); otherTp != other->trackPoints_.end();
         ++otherTp) {
      points.push_back(new TrackPoint(**otherTp, otherRepThisRep, &otherRepsToRemove));
    }
 
    insertPoints(points, id);
  }
 
 
  void Track::addTrackRep(AbsTrackRep* trackRep)
  {
    trackReps_.push_back(trackRep);
    fitStatuses_[trackRep] = new FitStatus();
  }
 
 
  void Track::deleteTrackRep(int id)
  {
    if (id < 0)
      id += trackReps_.size();
 
    AbsTrackRep* rep = trackReps_.at(id);
 
    // update cardinalRep_
    if (int(cardinalRep_) == id)
      cardinalRep_ = 0; // reset
    else if (int(cardinalRep_) > id)
      --cardinalRep_; // make cardinalRep_ point to the same TrackRep before and after deletion
 
    // delete FitterInfos related to the deleted TrackRep
    for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin(); pointIt != trackPoints_.end(); ++pointIt) {
      (*pointIt)->deleteFitterInfo(rep);
    }
 
    // delete fitStatus
    delete fitStatuses_.at(rep);
    fitStatuses_.erase(rep);
 
    // delete rep
    delete rep;
    trackReps_.erase(trackReps_.begin() + id);
  }
 
 
  void Track::setCardinalRep(int id)
  {
 
    if (id < 0)
      id += trackReps_.size();
 
    if (id >= 0 && (unsigned int)id < trackReps_.size())
      cardinalRep_ = id;
    else {
      cardinalRep_ = 0;
      errorOut << "Track::setCardinalRep: Attempted to set cardinalRep_ to a value out of bounds. Resetting  cardinalRep_ to 0." <<
               std::endl;
    }
  }
 
 
  void Track::determineCardinalRep()
  {
 
    // Todo: test
 
    if (trackReps_.size() <= 1)
      return;
 
    double minChi2(9.E99);
    const AbsTrackRep* bestRep(nullptr);
 
    for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it = fitStatuses_.begin(); it != fitStatuses_.end(); ++it) {
      if (it->second->isFitConverged()) {
        if (it->second->getChi2() < minChi2) {
          minChi2 = it->second->getChi2();
          bestRep = it->first;
        }
      }
    }
 
    if (bestRep != nullptr) {
      setCardinalRep(getIdForRep(bestRep));
    }
  }
 
 
  bool Track::sort()
  {
#ifdef DEBUG
    debugOut << "Track::sort \n";
#endif
 
    int nPoints(trackPoints_.size());
    // original order
    std::vector<TrackPoint*> pointsBefore(trackPoints_);
 
    // sort
    std::stable_sort(trackPoints_.begin(), trackPoints_.end(), TrackPointComparator());
 
    // see where order changed
    int equalUntil(-1), equalFrom(nPoints);
    for (int i = 0; i < nPoints; ++i) {
      if (pointsBefore[i] == trackPoints_[i])
        equalUntil = i;
      else
        break;
    }
 
    if (equalUntil == nPoints - 1)
      return false; // sorting did not change anything
 
 
    trackHasChanged();
 
    for (int i = nPoints - 1; i > equalUntil; --i) {
      if (pointsBefore[i] == trackPoints_[i])
        equalFrom = i;
      else
        break;
    }
 
#ifdef DEBUG
    debugOut << "Track::sort. Equal up to (including) hit " << equalUntil << " and from (including) hit " << equalFrom << " \n";
#endif
 
    deleteForwardInfo(equalUntil + 1, -1);
    deleteBackwardInfo(0, equalFrom - 1);
 
    deleteReferenceInfo(std::max(0, equalUntil + 1), std::min((int)trackPoints_.size() - 1, equalFrom - 1));
 
    fillPointsWithMeasurement();
 
    return true;
  }
 
 
  bool Track::udpateSeed(int id, AbsTrackRep* rep, bool biased)
  {
    try {
      const MeasuredStateOnPlane& fittedState = getFittedState(id, rep, biased);
      setTimeSeed(fittedState.getTime());
      setStateSeed(fittedState.get6DState());
      setCovSeed(fittedState.get6DCov());
 
      double fittedCharge = fittedState.getCharge();
 
      for (unsigned int i = 0; i < trackReps_.size(); ++i) {
        if (trackReps_[i]->getPDGCharge() * fittedCharge < 0) {
          trackReps_[i]->switchPDGSign();
        }
      }
    } catch (Exception& e) {
      // in this case the original track seed will be used
      return false;
    }
    return true;
  }
 
 
  void Track::reverseTrackPoints()
  {
 
    std::reverse(trackPoints_.begin(), trackPoints_.end());
 
    deleteForwardInfo(0, -1);
    deleteBackwardInfo(0, -1);
    deleteReferenceInfo(0, -1);
 
    fillPointsWithMeasurement();
  }
 
 
  void Track::switchPDGSigns(AbsTrackRep* rep)
  {
    if (rep != nullptr) {
      rep->switchPDGSign();
      return;
    }
 
    for (unsigned int i = 0; i < trackReps_.size(); ++i) {
      trackReps_[i]->switchPDGSign();
    }
  }
 
 
  void Track::reverseTrack()
  {
    udpateSeed(-1); // set fitted state of last hit as new seed
    reverseMomSeed(); // flip momentum direction
    switchPDGSigns();
    reverseTrackPoints(); // also deletes all fitterInfos
  }
 
 
  void Track::deleteForwardInfo(int startId, int endId, const AbsTrackRep* rep)
  {
#ifdef DEBUG
    debugOut << "Track::deleteForwardInfo from position " << startId  << " to " << endId << "\n";
#endif
 
    trackHasChanged();
 
    if (startId < 0)
      startId += trackPoints_.size();
    if (endId < 0)
      endId += trackPoints_.size();
    endId += 1;
 
    assert(endId >= startId);
 
    for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId;
         ++pointIt) {
      if (rep != nullptr) {
        if ((*pointIt)->hasFitterInfo(rep))
          (*pointIt)->getFitterInfo(rep)->deleteForwardInfo();
      } else {
        const std::vector<AbsFitterInfo*> fitterInfos = (*pointIt)->getFitterInfos();
        for (std::vector<AbsFitterInfo*>::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end();
             ++fitterInfoIt) {
          (*fitterInfoIt)->deleteForwardInfo();
        }
      }
    }
  }
 
  void Track::deleteBackwardInfo(int startId, int endId, const AbsTrackRep* rep)
  {
 
#ifdef DEBUG
    debugOut << "Track::deleteBackwardInfo from position " << startId  << " to " << endId << "\n";
#endif
 
    trackHasChanged();
 
    if (startId < 0)
      startId += trackPoints_.size();
    if (endId < 0)
      endId += trackPoints_.size();
    endId += 1;
 
    assert(endId >= startId);
 
 
    for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId;
         ++pointIt) {
      if (rep != nullptr) {
        if ((*pointIt)->hasFitterInfo(rep))
          (*pointIt)->getFitterInfo(rep)->deleteBackwardInfo();
      } else {
        const std::vector<AbsFitterInfo*> fitterInfos = (*pointIt)->getFitterInfos();
        for (std::vector<AbsFitterInfo*>::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end();
             ++fitterInfoIt) {
          (*fitterInfoIt)->deleteBackwardInfo();
        }
      }
    }
  }
 
  void Track::deleteReferenceInfo(int startId, int endId, const AbsTrackRep* rep)
  {
 
#ifdef DEBUG
    debugOut << "Track::deleteReferenceInfo from position " << startId  << " to " << endId << "\n";
#endif
 
    trackHasChanged();
 
    if (startId < 0)
      startId += trackPoints_.size();
    if (endId < 0)
      endId += trackPoints_.size();
    endId += 1;
 
    assert(endId >= startId);
 
    for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId;
         ++pointIt) {
      if (rep != nullptr) {
        if ((*pointIt)->hasFitterInfo(rep))
          (*pointIt)->getFitterInfo(rep)->deleteReferenceInfo();
      } else {
        std::vector<AbsFitterInfo*> fitterInfos = (*pointIt)->getFitterInfos();
        for (std::vector<AbsFitterInfo*>::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end();
             ++fitterInfoIt) {
          (*fitterInfoIt)->deleteReferenceInfo();
        }
      }
    }
  }
 
  void Track::deleteMeasurementInfo(int startId, int endId, const AbsTrackRep* rep)
  {
 
#ifdef DEBUG
    debugOut << "Track::deleteMeasurementInfo from position " << startId  << " to " << endId << "\n";
#endif
 
    trackHasChanged();
 
    if (startId < 0)
      startId += trackPoints_.size();
    if (endId < 0)
      endId += trackPoints_.size();
    endId += 1;
 
    assert(endId >= startId);
 
    for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId;
         ++pointIt) {
      if (rep != nullptr) {
        if ((*pointIt)->hasFitterInfo(rep))
          (*pointIt)->getFitterInfo(rep)->deleteMeasurementInfo();
      } else {
        std::vector<AbsFitterInfo*> fitterInfos = (*pointIt)->getFitterInfos();
        for (std::vector<AbsFitterInfo*>::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end();
             ++fitterInfoIt) {
          (*fitterInfoIt)->deleteMeasurementInfo();
        }
      }
    }
  }
 
  void Track::deleteFitterInfo(int startId, int endId, const AbsTrackRep* rep)
  {
 
#ifdef DEBUG
    debugOut << "Track::deleteFitterInfo from position " << startId  << " to " << endId << "\n";
#endif
 
    trackHasChanged();
 
    if (startId < 0)
      startId += trackPoints_.size();
    if (endId < 0)
      endId += trackPoints_.size();
    endId += 1;
 
    assert(endId >= startId);
 
    for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId;
         ++pointIt) {
      if (rep != nullptr) {
        if ((*pointIt)->hasFitterInfo(rep))
          (*pointIt)->deleteFitterInfo(rep);
      } else {
        for (std::vector<AbsTrackRep*>::const_iterator repIt = trackReps_.begin(); repIt != trackReps_.end(); ++repIt) {
          if ((*pointIt)->hasFitterInfo(*repIt))
            (*pointIt)->deleteFitterInfo(*repIt);
        }
      }
    }
  }
 
 
  double Track::getTrackLen(AbsTrackRep* rep, int startId, int endId) const
  {
 
    if (startId < 0)
      startId += trackPoints_.size();
    if (endId < 0)
      endId += trackPoints_.size();
 
    bool backwards(false);
    if (startId > endId) {
      double temp = startId;
      startId = endId;
      endId = temp;
      backwards = true;
    }
 
    endId += 1;
 
    if (rep == nullptr)
      rep = getCardinalRep();
 
    double trackLen(0);
    StateOnPlane state;
 
    for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId;
         ++pointIt) {
      if (!(*pointIt)->hasFitterInfo(rep)) {
        Exception e("Track::getTracklength: trackPoint has no fitterInfo", __LINE__, __FILE__);
        throw e;
      }
 
      if (pointIt != trackPoints_.begin() + startId) {
        trackLen += rep->extrapolateToPlane(state, (*pointIt)->getFitterInfo(rep)->getPlane());
      }
 
      state = (*pointIt)->getFitterInfo(rep)->getFittedState();
    }
 
    if (backwards)
      trackLen *= -1.;
 
    return trackLen;
  }
 
 
  TrackCand* Track::constructTrackCand() const
  {
    TrackCand* cand = new TrackCand();
 
    cand->setTime6DSeedAndPdgCode(timeSeed_, stateSeed_, getCardinalRep()->getPDG());
    cand->setCovSeed(covSeed_);
    cand->setMcTrackId(mcTrackId_);
 
    for (unsigned int i = 0; i < trackPointsWithMeasurement_.size(); ++i) {
      const TrackPoint* tp = trackPointsWithMeasurement_[i];
      const std::vector< AbsMeasurement* >& measurements = tp->getRawMeasurements();
 
      for (unsigned int j = 0; j < measurements.size(); ++j) {
        const AbsMeasurement* m = measurements[j];
        TrackCandHit* tch;
 
        int planeId = -1;
        if (dynamic_cast<const PlanarMeasurement*>(m)) {
          planeId = static_cast<const PlanarMeasurement*>(m)->getPlaneId();
        }
 
        if (m->isLeftRightMeasurement()) {
          tch = new WireTrackCandHit(m->getDetId(), m->getHitId(), planeId,
                                     tp->getSortingParameter(), m->getLeftRightResolution());
        } else {
          tch = new TrackCandHit(m->getDetId(), m->getHitId(), planeId,
                                 tp->getSortingParameter());
        }
        cand->addHit(tch);
      }
    }
 
    return cand;
  }
 
 
  double Track::getTOF(AbsTrackRep* rep, int startId, int endId) const
  {
 
    if (startId < 0)
      startId += trackPoints_.size();
    if (endId < 0)
      endId += trackPoints_.size();
 
    if (startId > endId) {
      std::swap(startId, endId);
    }
 
    endId += 1;
 
    if (rep == nullptr)
      rep = getCardinalRep();
 
    StateOnPlane state;
 
    const TrackPoint* startPoint(trackPoints_[startId]);
    const TrackPoint* endPoint(trackPoints_[endId]);
 
    if (!startPoint->hasFitterInfo(rep)
        || !endPoint->hasFitterInfo(rep)) {
      Exception e("Track::getTOF: trackPoint has no fitterInfo", __LINE__, __FILE__);
      throw e;
    }
 
    double tof = (rep->getTime(endPoint->getFitterInfo(rep)->getFittedState())
                  - rep->getTime(startPoint->getFitterInfo(rep)->getFittedState()));
    return tof;
  }
 
 
  void Track::fixWeights(AbsTrackRep* rep, int startId, int endId)
  {
 
    if (startId < 0)
      startId += trackPoints_.size();
    if (endId < 0)
      endId += trackPoints_.size();
 
    assert(startId >= 0);
    assert(startId <= endId);
    assert(endId <= (int)trackPoints_.size());
 
    std::vector< AbsFitterInfo* > fis;
 
    for (std::vector<TrackPoint*>::iterator tp = trackPoints_.begin() + startId; tp != trackPoints_.begin() + endId; ++tp) {
      fis.clear();
      if (rep == nullptr) {
        fis = (*tp)->getFitterInfos();
      } else if ((*tp)->hasFitterInfo(rep)) {
        fis.push_back((*tp)->getFitterInfo(rep));
      }
 
      for (std::vector< AbsFitterInfo* >::iterator fi = fis.begin(); fi != fis.end(); ++fi) {
        KalmanFitterInfo* kfi = dynamic_cast<KalmanFitterInfo*>(*fi);
        if (kfi == nullptr)
          continue;
 
        kfi->fixWeights();
      }
    }
  }
 
 
  void Track::prune(const Option_t* option)
  {
 
    PruneFlags f;
    f.setFlags(option);
 
    for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it = fitStatuses_.begin(); it != fitStatuses_.end(); ++it) {
      it->second->getPruneFlags().setFlags(option);
    }
 
    // prune trackPoints
    if (f.hasFlags("F") || f.hasFlags("L")) {
      TrackPoint* firstPoint = getPointWithFitterInfo(0);
      TrackPoint* lastPoint = getPointWithFitterInfo(-1);
      for (unsigned int i = 0; i < trackPoints_.size(); ++i) {
        if (trackPoints_[i] == firstPoint && f.hasFlags("F"))
          continue;
 
        if (trackPoints_[i] == lastPoint && f.hasFlags("L"))
          continue;
 
        delete trackPoints_[i];
        trackPoints_.erase(trackPoints_.begin() + i);
        --i;
      }
    }
 
    // prune TrackReps
    if (f.hasFlags("C")) {
      for (unsigned int i = 0; i < trackReps_.size(); ++i) {
        if (i != cardinalRep_) {
          deleteTrackRep(i);
          --i;
        }
      }
    }
 
 
    // from remaining trackPoints: prune measurementsOnPlane, unneeded fitterInfoStuff
    for (unsigned int i = 0; i < trackPoints_.size(); ++i) {
      if (f.hasFlags("W"))
        trackPoints_[i]->deleteRawMeasurements();
 
      std::vector< AbsFitterInfo* > fis =  trackPoints_[i]->getFitterInfos();
      for (unsigned int j = 0; j < fis.size(); ++j) {
 
        if (i == 0 && f.hasFlags("FLI"))
          fis[j]->deleteForwardInfo();
        else if (i == trackPoints_.size() - 1 && f.hasFlags("FLI"))
          fis[j]->deleteBackwardInfo();
        else if (f.hasFlags("FI"))
          fis[j]->deleteForwardInfo();
        else if (f.hasFlags("LI"))
          fis[j]->deleteBackwardInfo();
 
        if (f.hasFlags("U") && dynamic_cast<KalmanFitterInfo*>(fis[j]) != nullptr) {
          static_cast<KalmanFitterInfo*>(fis[j])->deletePredictions();
        }
 
        // also delete reference info if points have been removed since it is invalid then!
        if (f.hasFlags("R") or f.hasFlags("F") or f.hasFlags("L"))
          fis[j]->deleteReferenceInfo();
        if (f.hasFlags("M"))
          fis[j]->deleteMeasurementInfo();
      }
    }
 
    fillPointsWithMeasurement();
 
#ifdef DEBUG
    debugOut << "pruned Track: "; Print();
#endif
 
  }
 
 
  void Track::Print(const Option_t* option) const
  {
    TString opt = option;
    opt.ToUpper();
    if (opt.Contains("C")) { // compact
 
      printOut << "\n    ";
      for (unsigned int i = 0; i < trackPoints_.size(); ++i) {
 
        int color = 32 * (size_t)(trackPoints_[i]) % 15;
        switch (color) {
          case 0:
            printOut << "\033[1;30m";
            break;
          case 1:
            printOut << "\033[0;34m";
            break;
          case 2:
            printOut << "\033[1;34m";
            break;
          case 3:
            printOut << "\033[0;32m";
            break;
          case 4:
            printOut << "\033[1;32m";
            break;
          case 5:
            printOut << "\033[0;36m";
            break;
          case 6:
            printOut << "\033[1;36m";
            break;
          case 7:
            printOut << "\033[0;31m";
            break;
          case 8:
            printOut << "\033[1;31m";
            break;
          case 9:
            printOut << "\033[0;35m";
            break;
          case 10:
            printOut << "\033[1;35m";
            break;
          case 11:
            printOut << "\033[0;33m";
            break;
          case 12:
            printOut << "\033[1;33m";
            break;
          case 13:
            printOut << "\033[0;37m";
            break;
          default:
            ;
        }
        printOut << trackPoints_[i] << "\033[00m  ";
      }
      printOut << "\n";
 
      printOut << "   ";
      for (unsigned int i = 0; i < trackPoints_.size(); ++i) {
        printf("% -9.3g  ", trackPoints_[i]->getSortingParameter());
      }
 
      for (std::vector<AbsTrackRep*>::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) {
        printOut << "\n" << getIdForRep(*rep) << "   ";
        for (unsigned int i = 0; i < trackPoints_.size(); ++i) {
          if (! trackPoints_[i]->hasFitterInfo(*rep)) {
            printOut << "           ";
            continue;
          }
          AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep);
          if (fi->hasMeasurements())
            printOut << "M";
          else
            printOut << " ";
 
          if (fi->hasReferenceState())
            printOut << "R";
          else
            printOut << " ";
 
          printOut << "         ";
        }
        printOut << "\n";
 
        printOut << " -> ";
        for (unsigned int i = 0; i < trackPoints_.size(); ++i) {
          if (! trackPoints_[i]->hasFitterInfo(*rep)) {
            printOut << "           ";
            continue;
          }
          AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep);
          if (fi->hasForwardPrediction())
            printOut << "P";
          else
            printOut << " ";
 
          if (fi->hasForwardUpdate())
            printOut << "U";
          else
            printOut << " ";
 
          printOut << "         ";
        }
        printOut << "\n";
 
        printOut << " <- ";
        for (unsigned int i = 0; i < trackPoints_.size(); ++i) {
          if (! trackPoints_[i]->hasFitterInfo(*rep)) {
            printOut << "           ";
            continue;
          }
          AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep);
          if (fi->hasBackwardPrediction())
            printOut << "P";
          else
            printOut << " ";
 
          if (fi->hasBackwardUpdate())
            printOut << "U";
          else
            printOut << " ";
 
          printOut << "         ";
        }
 
        printOut << "\n";
 
      } //end loop over reps
 
      printOut << "\n";
      return;
    }
 
 
 
    printOut << "=======================================================================================\n";
    printOut << "genfit::Track, containing " << trackPoints_.size() << " TrackPoints and " << trackReps_.size() << " TrackReps.\n";
    printOut << " Seed state: "; stateSeed_.Print();
 
    for (unsigned int i = 0; i < trackReps_.size(); ++i) {
      printOut << " TrackRep Nr. " << i;
      if (i == cardinalRep_)
        printOut << " (This is the cardinal rep)";
      printOut << "\n";
      trackReps_[i]->Print();
    }
 
    printOut << "---------------------------------------------------------------------------------------\n";
 
    for (unsigned int i = 0; i < trackPoints_.size(); ++i) {
      printOut << "TrackPoint Nr. " << i << "\n";
      trackPoints_[i]->Print();
      printOut << "..........................................................................\n";
    }
 
    for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it = fitStatuses_.begin(); it != fitStatuses_.end(); ++it) {
      it->second->Print();
    }
 
    printOut << "=======================================================================================\n";
 
  }
 
 
  void Track::checkConsistency() const
  {
 
    // cppcheck-suppress unreadVariable
    bool consistent = true;
    std::stringstream failures;
 
    std::map<const AbsTrackRep*, const KalmanFitterInfo*> prevFis;
 
    // check if seed is 6D
    if (stateSeed_.GetNrows() != 6) {
      failures << "Track::checkConsistency(): stateSeed_ dimension != 6" << std::endl;
      // cppcheck-suppress unreadVariable
      consistent = false;
    }
 
    if (covSeed_.GetNrows() != 6) {
      failures << "Track::checkConsistency(): covSeed_ dimension != 6" << std::endl;
      // cppcheck-suppress unreadVariable
      consistent = false;
    }
 
    if (covSeed_.Max() == 0.) {
      // Nota bene: The consistency is not set to false when this occurs, because it does not break the consistency of
      // the track. However, when something else fails we keep this as additional error information.
      failures << "Track::checkConsistency(): Warning: covSeed_ is zero" << std::endl;
    }
 
    // check if correct number of fitStatuses
    if (fitStatuses_.size() != trackReps_.size()) {
      failures << "Track::checkConsistency(): Number of fitStatuses is != number of TrackReps " << std::endl;
      // cppcheck-suppress unreadVariable
      consistent = false;
    }
 
    // check if cardinalRep_ is in range of trackReps_
    if (trackReps_.size() && cardinalRep_ >= trackReps_.size()) {
      failures << "Track::checkConsistency(): cardinalRep id " << cardinalRep_ << " out of bounds" << std::endl;
      // cppcheck-suppress unreadVariable
      consistent = false;
    }
 
    for (std::vector<AbsTrackRep*>::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) {
      // check for nullptr
      if ((*rep) == nullptr) {
        failures << "Track::checkConsistency(): TrackRep is nullptr" << std::endl;
        // cppcheck-suppress unreadVariable
        consistent = false;
      }
 
      // check for valid pdg code
      TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle((*rep)->getPDG());
      if (particle == nullptr) {
        failures << "Track::checkConsistency(): TrackRep pdg ID " << (*rep)->getPDG() << " is not valid" << std::endl;
        // cppcheck-suppress unreadVariable
        consistent = false;
      }
 
      // check if corresponding FitStatus is there
      if (fitStatuses_.find(*rep) == fitStatuses_.end() and fitStatuses_.find(*rep)->second != nullptr) {
        failures << "Track::checkConsistency(): No FitStatus for Rep or FitStatus is nullptr" << std::endl;
        // cppcheck-suppress unreadVariable
        consistent = false;
      }
    }
 
    // check TrackPoints
    for (std::vector<TrackPoint*>::const_iterator tp = trackPoints_.begin(); tp != trackPoints_.end(); ++tp) {
      // check for nullptr
      if ((*tp) == nullptr) {
        failures << "Track::checkConsistency(): TrackPoint is nullptr" << std::endl;
        // cppcheck-suppress unreadVariable
        consistent = false;
      }
      // check if trackPoint points back to this track
      if ((*tp)->getTrack() != this) {
        failures << "Track::checkConsistency(): TrackPoint does not point back to this track" << std::endl;
        // cppcheck-suppress unreadVariable
        consistent = false;
      }
 
      // check rawMeasurements
      const std::vector<AbsMeasurement*>& rawMeasurements = (*tp)->getRawMeasurements();
      for (std::vector<AbsMeasurement*>::const_iterator m = rawMeasurements.begin(); m != rawMeasurements.end(); ++m) {
        // check for nullptr
        if ((*m) == nullptr) {
          failures << "Track::checkConsistency(): Measurement is nullptr" << std::endl;
          // cppcheck-suppress unreadVariable
          consistent = false;
        }
        // check if measurement points back to TrackPoint
        if ((*m)->getTrackPoint() != *tp) {
          failures << "Track::checkConsistency(): Measurement does not point back to correct TrackPoint" << std::endl;
          // cppcheck-suppress unreadVariable
          consistent = false;
        }
      }
 
      // check fitterInfos
      std::vector<AbsFitterInfo*> fitterInfos = (*tp)->getFitterInfos();
      for (std::vector<AbsFitterInfo*>::const_iterator fi = fitterInfos.begin(); fi != fitterInfos.end(); ++fi) {
        // check for nullptr
        if ((*fi) == nullptr) {
          failures << "Track::checkConsistency(): FitterInfo is nullptr. TrackPoint: " << *tp << std::endl;
          // cppcheck-suppress unreadVariable
          consistent = false;
        }
 
        // check if fitterInfos point to valid TrackReps in trackReps_
        int mycount(0);
        for (std::vector<AbsTrackRep*>::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) {
          if ((*rep) == (*fi)->getRep()) {
            ++mycount;
          }
        }
        if (mycount ==  0) {
          failures << "Track::checkConsistency(): fitterInfo points to TrackRep which is not in Track" << std::endl;
          // cppcheck-suppress unreadVariable
          consistent = false;
        }
 
        if (!((*fi)->checkConsistency(&(this->getFitStatus((*fi)->getRep())->getPruneFlags())))) {
          failures << "Track::checkConsistency(): FitterInfo not consistent. TrackPoint: " << *tp << std::endl;
          // cppcheck-suppress unreadVariable
          consistent = false;
        }
 
        if (dynamic_cast<KalmanFitterInfo*>(*fi) != nullptr) {
          if (prevFis[(*fi)->getRep()] != nullptr &&
              static_cast<KalmanFitterInfo*>(*fi)->hasReferenceState() &&
              prevFis[(*fi)->getRep()]->hasReferenceState()) {
            double len = static_cast<KalmanFitterInfo*>(*fi)->getReferenceState()->getForwardSegmentLength();
            double prevLen = prevFis[(*fi)->getRep()]->getReferenceState()->getBackwardSegmentLength();
            if (fabs(prevLen + len) > 1E-10) {
              failures << "Track::checkConsistency(): segment lengths of reference states for rep " << (*fi)->getRep() << " (id " << getIdForRep((
                         *fi)->getRep()) << ") at TrackPoint " << (*tp) << " don't match" << std::endl;
              failures << prevLen << " + " << len << " = " << prevLen + len << std::endl;
              failures << "TrackPoint " << *tp << ", FitterInfo " << *fi << ", rep " << getIdForRep((*fi)->getRep()) << std::endl;
              // cppcheck-suppress unreadVariable
              consistent = false;
            }
          }
 
          prevFis[(*fi)->getRep()] = static_cast<KalmanFitterInfo*>(*fi);
        } else
          prevFis[(*fi)->getRep()] = nullptr;
 
      } // end loop over FitterInfos
 
    } // end loop over TrackPoints
 
 
    // check trackPointsWithMeasurement_
    std::vector<TrackPoint*> trackPointsWithMeasurement;
 
    for (std::vector<TrackPoint*>::const_iterator it = trackPoints_.begin(); it != trackPoints_.end(); ++it) {
      if ((*it)->hasRawMeasurements()) {
        trackPointsWithMeasurement.push_back(*it);
      }
    }
 
    if (trackPointsWithMeasurement.size() != trackPointsWithMeasurement_.size()) {
      failures << "Track::checkConsistency(): trackPointsWithMeasurement_ has incorrect size" << std::endl;
      // cppcheck-suppress unreadVariable
      consistent = false;
    }
 
    for (unsigned int i = 0; i < trackPointsWithMeasurement.size(); ++i) {
      if (trackPointsWithMeasurement[i] != trackPointsWithMeasurement_[i]) {
        failures << "Track::checkConsistency(): trackPointsWithMeasurement_ is not correct" << std::endl;
        failures << "has         id " << i << ", address " << trackPointsWithMeasurement_[i] << std::endl;
        failures << "should have id " << i << ", address " << trackPointsWithMeasurement[i] << std::endl;
        // cppcheck-suppress unreadVariable
        consistent = false;
      }
    }
 
    if (not consistent) {
      throw genfit::Exception(failures.str(), __LINE__, __FILE__);
    }
  }
 
 
  void Track::trackHasChanged()
  {
 
#ifdef DEBUG
    debugOut << "Track::trackHasChanged \n";
#endif
 
    if (fitStatuses_.empty())
      return;
 
    for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it = fitStatuses_.begin(); it != fitStatuses_.end(); ++it) {
      it->second->setHasTrackChanged();
    }
  }
 
 
  void Track::fillPointsWithMeasurement()
  {
    trackPointsWithMeasurement_.clear();
    trackPointsWithMeasurement_.reserve(trackPoints_.size());
 
    for (std::vector<TrackPoint*>::const_iterator it = trackPoints_.begin(); it != trackPoints_.end(); ++it) {
      if ((*it)->hasRawMeasurements()) {
        trackPointsWithMeasurement_.push_back(*it);
      }
    }
  }
 
 
  void Track::Streamer(TBuffer& R__b)
  {
    // Stream an object of class genfit::Track.
    const bool streamTrackPoints = true; // debugging option
    //This works around a msvc bug and should be harmless on other platforms
    typedef ::genfit::Track thisClass;
    UInt_t R__s, R__c;
    if (R__b.IsReading()) {
      this->Clear();
      Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
      if (R__v) { }
      TObject::Streamer(R__b);
      {
        std::vector<AbsTrackRep*>& R__stl =  trackReps_;
        R__stl.clear();
        TClass* R__tcl1 = TBuffer::GetClass(typeid(genfit::AbsTrackRep));
        if (R__tcl1 == 0) {
          Error("trackReps_ streamer", "Missing the TClass object for genfit::AbsTrackRep!");
          return;
        }
        int R__i, R__n;
        R__b >> R__n;
        R__stl.reserve(R__n);
        for (R__i = 0; R__i < R__n; R__i++) {
          genfit::AbsTrackRep* R__t;
          R__b >> R__t;
          R__stl.push_back(R__t);
        }
      }
      R__b >> cardinalRep_;
      if (streamTrackPoints) {
        std::vector<TrackPoint*>& R__stl =  trackPoints_;
        R__stl.clear();
        TClass* R__tcl1 = TBuffer::GetClass(typeid(genfit::TrackPoint));
        if (R__tcl1 == 0) {
          Error("trackPoints_ streamer", "Missing the TClass object for genfit::TrackPoint!");
          return;
        }
        int R__i, R__n;
        R__b >> R__n;
        R__stl.reserve(R__n);
        for (R__i = 0; R__i < R__n; R__i++) {
          genfit::TrackPoint* R__t;
          R__t = (genfit::TrackPoint*)R__b.ReadObjectAny(R__tcl1);
          R__t->setTrack(this);
          R__t->fixupRepsForReading();
          R__stl.push_back(R__t);
        }
      }
      {
        std::map<const AbsTrackRep*, FitStatus*>& R__stl =  fitStatuses_;
        R__stl.clear();
        TClass* R__tcl1 = TBuffer::GetClass(typeid(genfit::AbsTrackRep));
        if (R__tcl1 == 0) {
          Error("fitStatuses_ streamer", "Missing the TClass object for genfit::AbsTrackRep!");
          return;
        }
        TClass* R__tcl2 = TBuffer::GetClass(typeid(genfit::FitStatus));
        if (R__tcl2 == 0) {
          Error("fitStatuses_ streamer", "Missing the TClass object for genfit::FitStatus!");
          return;
        }
        int R__i, R__n;
        R__b >> R__n;
        for (R__i = 0; R__i < R__n; R__i++) {
          int id;
          R__b >> id;
          genfit::FitStatus* R__t2;
          R__t2 = (genfit::FitStatus*)R__b.ReadObjectAny(R__tcl2);
 
          R__stl[this->getTrackRep(id)] = R__t2;
        }
      }
      // timeSeed_ was introduced in version 3.  If reading an earlier
      // version, default to zero.
      if (R__v >= 3) {
        R__b >> timeSeed_;
      } else {
        timeSeed_ = 0;
      }
      stateSeed_.Streamer(R__b);
      covSeed_.Streamer(R__b);
      R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
 
      fillPointsWithMeasurement();
    } else {
      R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
      TObject::Streamer(R__b);
      {
        std::vector<AbsTrackRep*>& R__stl =  trackReps_;
        int R__n = int(R__stl.size());
        R__b << R__n;
        if (R__n) {
          TClass* R__tcl1 = TBuffer::GetClass(typeid(genfit::AbsTrackRep));
          if (R__tcl1 == 0) {
            Error("trackReps_ streamer", "Missing the TClass object for genfit::AbsTrackRep!");
            return;
          }
          std::vector<AbsTrackRep*>::iterator R__k;
          for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
            R__b << *R__k;
          }
        }
      }
      R__b << cardinalRep_;
      if (streamTrackPoints) {
        std::vector<TrackPoint*>& R__stl =  trackPoints_;
        int R__n = int(R__stl.size());
        R__b << R__n;
        if (R__n) {
          std::vector<TrackPoint*>::iterator R__k;
          for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
            R__b << (*R__k);
          }
        }
      }
      {
        std::map<const AbsTrackRep*, FitStatus*>& R__stl =  fitStatuses_;
        int R__n = int(R__stl.size());
        R__b << R__n;
        if (R__n) {
          TClass* R__tcl1 = TBuffer::GetClass(typeid(genfit::AbsTrackRep));
          if (R__tcl1 == 0) {
            Error("fitStatuses_ streamer", "Missing the TClass object for genfit::AbsTrackRep!");
            return;
          }
          TClass* R__tcl2 = TBuffer::GetClass(typeid(genfit::FitStatus));
          if (R__tcl2 == 0) {
            Error("fitStatuses_ streamer", "Missing the TClass object for genfit::FitStatus!");
            return;
          }
          std::map<const AbsTrackRep*, FitStatus*>::iterator R__k;
          for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
            int id = this->getIdForRep((*R__k).first);
            R__b << id;
            R__b << ((*R__k).second);
          }
        }
      }
      R__b << timeSeed_;
      stateSeed_.Streamer(R__b);
      covSeed_.Streamer(R__b);
      R__b.SetByteCount(R__c, kTRUE);
    }
  }
 
  void Track::deleteTrackPointsAndFitStatus()
  {
    for (size_t i = 0; i < trackPoints_.size(); ++i)
      delete trackPoints_[i];
 
    trackPoints_.clear();
    trackPointsWithMeasurement_.clear();
 
    for (std::map< const AbsTrackRep*, FitStatus* >::iterator it = fitStatuses_.begin(); it != fitStatuses_.end(); ++it)
      delete it->second;
    fitStatuses_.clear();
  }
 
} /* End of namespace genfit */