/* 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 . */ #include "Track.h" #include "Exception.h" #include "KalmanFitterInfo.h" #include "KalmanFitStatus.h" #include "PlanarMeasurement.h" #include "AbsMeasurement.h" #include "WireTrackCandHit.h" #include #include #include #include #include #include "TBuffer.h" //#include //#define DEBUG namespace genfit { Track::Track() : cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(6), covSeed_(6) { ; } Track::Track(const TrackCand& trackCand, const MeasurementFactory& factory, AbsTrackRep* rep) : cardinalRep_(0), fitStatuses_(), stateSeed_(6), covSeed_(6) { if (rep != NULL) addTrackRep(rep); // create the measurements using the factory. std::vector 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; isetSortingParameter(trackCand.getHit(i)->getSortingParameter()); insertPoint(tp); } // Copy seed information from candidate timeSeed_ = trackCand.getTimeSeed(); stateSeed_ = trackCand.getStateSeed(); covSeed_ = trackCand.getCovSeed(); mcTrackId_ = trackCand.getMcTrackId(); // fill cache fillPointsWithMeasurement(); // self test assert(checkConsistency()); } 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_) { assert(rhs.checkConsistency()); std::map oldRepNewRep; for (std::vector::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::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(); // self test assert(checkConsistency()); } Track& Track::operator=(Track other) { swap(other); for (std::vector::const_iterator it=trackPoints_.begin(); it!=trackPoints_.end(); ++it) { trackPoints_.back()->setTrack(this); } fillPointsWithMeasurement(); // self test assert(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 { int i(0); for (std::vector::const_iterator it = trackPointsWithMeasurement_.begin(); it != trackPointsWithMeasurement_.end(); ++it) { if ((*it)->hasFitterInfo(rep)) { if (id == i) return (*it); ++i; } } if (i == 0) return NULL; if (id < 0) id += i; return getPointWithMeasurementAndFitterInfo(id, rep); } TrackPoint* Track::getPointWithFitterInfo(int id, const AbsTrackRep* rep) const { if (rep == NULL) rep = getCardinalRep(); if (id >= 0) { int i = 0; for (std::vector::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::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 == NULL) rep = getCardinalRep(); TrackPoint* point = getPointWithFitterInfo(id, rep); if (point == NULL) { 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 == NULL) rep = getCardinalRep(); if (fitStatuses_.find(rep) == fitStatuses_.end()) return false; return (fitStatuses_.at(rep) != NULL); } bool Track::hasKalmanFitStatus(const AbsTrackRep* rep) const { if (rep == NULL) rep = getCardinalRep(); if (fitStatuses_.find(rep) == fitStatuses_.end()) return false; return (dynamic_cast(fitStatuses_.at(rep)) != NULL); } KalmanFitStatus* Track::getKalmanFitStatus(const AbsTrackRep* rep) const { return dynamic_cast(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 std::cout << "Track::insertPoint at position " << id << "\n"; #endif assert(point!=NULL); 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 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::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 std::cout << "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::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::mergeTrack(const Track* other, int id) { #ifdef DEBUG std::cout << "Track::mergeTrack\n"; #endif if (other->getNumPoints() == 0) return; std::map otherRepThisRep; std::vector otherRepsToRemove; for (std::vector::const_iterator otherRep=other->trackReps_.begin(); otherRep!=other->trackReps_.end(); ++otherRep) { bool found(false); for (std::vector::const_iterator thisRep=trackReps_.begin(); thisRep!=trackReps_.end(); ++thisRep) { if ((*thisRep)->isSame(*otherRep)) { otherRepThisRep[*otherRep] = *thisRep; #ifdef DEBUG std::cout << " 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 std::cout << " remove other rep " << *otherRep << "\n"; #endif } } std::vector points; points.reserve(other->getNumPoints()); for (std::vector::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::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; std::cerr << "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(NULL); 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 != NULL) { setCardinalRep(getIdForRep(bestRep)); } } bool Track::sort() { #ifdef DEBUG std::cout << "Track::sort \n"; #endif int nPoints(trackPoints_.size()); // original order std::vector pointsBefore(trackPoints_); // sort std::stable_sort(trackPoints_.begin(), trackPoints_.end(), TrackPointComparator()); // see where order changed int equalUntil(-1), equalFrom(nPoints); for (int i = 0; iequalUntil; --i) { if (pointsBefore[i] == trackPoints_[i]) equalFrom = i; else break; } #ifdef DEBUG std::cout << "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; igetPDGCharge() * 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 != NULL) { rep->switchPDGSign(); return; } for (unsigned int i = 0; iswitchPDGSign(); } } 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 std::cout << "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::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) { if (rep != NULL) { if ((*pointIt)->hasFitterInfo(rep)) (*pointIt)->getFitterInfo(rep)->deleteForwardInfo(); } else { const std::vector fitterInfos = (*pointIt)->getFitterInfos(); for (std::vector::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end(); ++fitterInfoIt) { (*fitterInfoIt)->deleteForwardInfo(); } } } } void Track::deleteBackwardInfo(int startId, int endId, const AbsTrackRep* rep) { #ifdef DEBUG std::cout << "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::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) { if (rep != NULL) { if ((*pointIt)->hasFitterInfo(rep)) (*pointIt)->getFitterInfo(rep)->deleteBackwardInfo(); } else { const std::vector fitterInfos = (*pointIt)->getFitterInfos(); for (std::vector::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end(); ++fitterInfoIt) { (*fitterInfoIt)->deleteBackwardInfo(); } } } } void Track::deleteReferenceInfo(int startId, int endId, const AbsTrackRep* rep) { #ifdef DEBUG std::cout << "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::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) { if (rep != NULL) { if ((*pointIt)->hasFitterInfo(rep)) (*pointIt)->getFitterInfo(rep)->deleteReferenceInfo(); } else { std::vector fitterInfos = (*pointIt)->getFitterInfos(); for (std::vector::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end(); ++fitterInfoIt) { (*fitterInfoIt)->deleteReferenceInfo(); } } } } void Track::deleteMeasurementInfo(int startId, int endId, const AbsTrackRep* rep) { #ifdef DEBUG std::cout << "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::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) { if (rep != NULL) { if ((*pointIt)->hasFitterInfo(rep)) (*pointIt)->getFitterInfo(rep)->deleteMeasurementInfo(); } else { std::vector fitterInfos = (*pointIt)->getFitterInfos(); for (std::vector::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end(); ++fitterInfoIt) { (*fitterInfoIt)->deleteMeasurementInfo(); } } } } void Track::deleteFitterInfo(int startId, int endId, const AbsTrackRep* rep) { #ifdef DEBUG std::cout << "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::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) { if (rep != NULL) { if ((*pointIt)->hasFitterInfo(rep)) (*pointIt)->deleteFitterInfo(rep); } else { for (std::vector::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 == NULL) rep = getCardinalRep(); double trackLen(0); StateOnPlane state; for (std::vector::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(m)) { planeId = static_cast(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 == NULL) 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::iterator tp = trackPoints_.begin() + startId; tp != trackPoints_.begin() + endId; ++tp) { fis.clear(); if (rep == NULL) { 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(*fi); if (kfi == NULL) 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; ideleteRawMeasurements(); std::vector< AbsFitterInfo* > fis = trackPoints_[i]->getFitterInfos(); for (unsigned int j = 0; jdeleteForwardInfo(); 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(fis[j]) != NULL) { static_cast(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 std::cout << "pruned Track: "; Print(); #endif } void Track::Print(const Option_t* option) const { TString opt = option; opt.ToUpper(); if (opt.Contains("C")) { // compact std::cout << "\n "; for (unsigned int i=0; igetSortingParameter()); } for (std::vector::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) { std::cout << "\n" << getIdForRep(*rep) << " "; for (unsigned int i=0; ihasFitterInfo(*rep)) { std::cout << " "; continue; } AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep); if (fi->hasMeasurements()) std::cout << "M"; else std::cout << " "; if (fi->hasReferenceState()) std::cout << "R"; else std::cout << " "; std::cout << " "; } std::cout << "\n"; std::cout << " -> "; for (unsigned int i=0; ihasFitterInfo(*rep)) { std::cout << " "; continue; } AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep); if (fi->hasForwardPrediction()) std::cout << "P"; else std::cout << " "; if (fi->hasForwardUpdate()) std::cout << "U"; else std::cout << " "; std::cout << " "; } std::cout << "\n"; std::cout << " <- "; for (unsigned int i=0; ihasFitterInfo(*rep)) { std::cout << " "; continue; } AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep); if (fi->hasBackwardPrediction()) std::cout << "P"; else std::cout << " "; if (fi->hasBackwardUpdate()) std::cout << "U"; else std::cout << " "; std::cout << " "; } std::cout << "\n"; } //end loop over reps std::cout << "\n"; return; } std::cout << "=======================================================================================\n"; std::cout << "genfit::Track, containing " << trackPoints_.size() << " TrackPoints and " << trackReps_.size() << " TrackReps.\n"; std::cout << " Seed state: "; stateSeed_.Print(); for (unsigned int i=0; iPrint(); } std::cout << "---------------------------------------------------------------------------------------\n"; for (unsigned int i=0; iPrint(); std::cout << "..........................................................................\n"; } for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it=fitStatuses_.begin(); it!=fitStatuses_.end(); ++it) { it->second->Print(); } std::cout << "=======================================================================================\n"; } bool Track::checkConsistency() const { bool retVal(true); std::map prevFis; // check if seed is 6D if (stateSeed_.GetNrows() != 6) { std::cerr << "Track::checkConsistency(): stateSeed_ dimension != 6" << std::endl; retVal = false; } if (covSeed_.GetNrows() != 6) { std::cerr << "Track::checkConsistency(): covSeed_ dimension != 6" << std::endl; retVal = false; } if (covSeed_.Max() == 0.) { std::cerr << "Track::checkConsistency(): Warning: covSeed_ is zero" << std::endl; //retVal = false; } // check if correct number of fitStatuses if (fitStatuses_.size() != trackReps_.size()) { std::cerr << "Track::checkConsistency(): Number of fitStatuses is != number of TrackReps " << std::endl; retVal = false; } // check if cardinalRep_ is in range of trackReps_ if (trackReps_.size() && cardinalRep_ >= trackReps_.size()) { std::cerr << "Track::checkConsistency(): cardinalRep id " << cardinalRep_ << " out of bounds" << std::endl; retVal = false; } for (std::vector::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) { // check for NULL if ((*rep) == NULL) { std::cerr << "Track::checkConsistency(): TrackRep is NULL" << std::endl; retVal = false; } // check for valid pdg code TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle((*rep)->getPDG()); if (particle == NULL) { std::cerr << "Track::checkConsistency(): TrackRep pdg ID " << (*rep)->getPDG() << " is not valid" << std::endl; retVal = false; } // check if corresponding FitStatus is there if (fitStatuses_.find(*rep) == fitStatuses_.end() and fitStatuses_.find(*rep)->second != NULL) { std::cerr << "Track::checkConsistency(): No FitStatus for Rep or FitStatus is NULL" << std::endl; retVal = false; } } // check TrackPoints for (std::vector::const_iterator tp = trackPoints_.begin(); tp != trackPoints_.end(); ++tp) { // check for NULL if ((*tp) == NULL) { std::cerr << "Track::checkConsistency(): TrackPoint is NULL" << std::endl; retVal = false; } // check if trackPoint points back to this track if ((*tp)->getTrack() != this) { std::cerr << "Track::checkConsistency(): TrackPoint does not point back to this track" << std::endl; retVal = false; } // check rawMeasurements const std::vector& rawMeasurements = (*tp)->getRawMeasurements(); for (std::vector::const_iterator m = rawMeasurements.begin(); m != rawMeasurements.end(); ++m) { // check for NULL if ((*m) == NULL) { std::cerr << "Track::checkConsistency(): Measurement is NULL" << std::endl; retVal = false; } // check if measurement points back to TrackPoint if ((*m)->getTrackPoint() != *tp) { std::cerr << "Track::checkConsistency(): Measurement does not point back to correct TrackPoint" << std::endl; retVal = false; } } // check fitterInfos std::vector fitterInfos = (*tp)->getFitterInfos(); for (std::vector::const_iterator fi = fitterInfos.begin(); fi != fitterInfos.end(); ++fi) { // check for NULL if ((*fi) == NULL) { std::cerr << "Track::checkConsistency(): FitterInfo is NULL. TrackPoint: " << *tp << std::endl; retVal = false; } // check if fitterInfos point to valid TrackReps in trackReps_ int mycount (0); for (std::vector::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) { if ( (*rep) == (*fi)->getRep() ) { ++mycount; } } if (mycount == 0) { std::cerr << "Track::checkConsistency(): fitterInfo points to TrackRep which is not in Track" << std::endl; retVal = false; } if (!( (*fi)->checkConsistency(&(this->getFitStatus((*fi)->getRep())->getPruneFlags())) ) ) { std::cerr << "Track::checkConsistency(): FitterInfo not consistent. TrackPoint: " << *tp << std::endl; retVal = false; } if (dynamic_cast(*fi) != NULL) { if (prevFis[(*fi)->getRep()] != NULL && static_cast(*fi)->hasReferenceState() && prevFis[(*fi)->getRep()]->hasReferenceState() ) { double len = static_cast(*fi)->getReferenceState()->getForwardSegmentLength(); double prevLen = prevFis[(*fi)->getRep()]->getReferenceState()->getBackwardSegmentLength(); if (fabs(prevLen + len) > 1E-10 ) { std::cerr << "Track::checkConsistency(): segment lengths of reference states for rep " << (*fi)->getRep() << " (id " << getIdForRep((*fi)->getRep()) << ") at TrackPoint " << (*tp) << " don't match" << std::endl; std::cerr << prevLen << " + " << len << " = " << prevLen + len << std::endl; std::cerr << "TrackPoint " << *tp << ", FitterInfo " << *fi << ", rep " << getIdForRep((*fi)->getRep()) << std::endl; retVal = false; } } prevFis[(*fi)->getRep()] = static_cast(*fi); } else prevFis[(*fi)->getRep()] = NULL; } // end loop over FitterInfos } // end loop over TrackPoints // check trackPointsWithMeasurement_ std::vector trackPointsWithMeasurement; for (std::vector::const_iterator it = trackPoints_.begin(); it != trackPoints_.end(); ++it) { if ((*it)->hasRawMeasurements()) { trackPointsWithMeasurement.push_back(*it); } } if (trackPointsWithMeasurement.size() != trackPointsWithMeasurement_.size()) { std::cerr << "Track::checkConsistency(): trackPointsWithMeasurement_ has incorrect size" << std::endl; retVal = false; } for (unsigned int i = 0; i < trackPointsWithMeasurement.size(); ++i) { if (trackPointsWithMeasurement[i] != trackPointsWithMeasurement_[i]) { std::cerr << "Track::checkConsistency(): trackPointsWithMeasurement_ is not correct" << std::endl; std::cerr << "has id " << i << ", address " << trackPointsWithMeasurement_[i] << std::endl; std::cerr << "should have id " << i << ", address " << trackPointsWithMeasurement[i] << std::endl; retVal = false; } } return retVal; } void Track::trackHasChanged() { #ifdef DEBUG std::cout << "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::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 &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 &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 &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 &R__stl = trackReps_; int R__n=(&R__stl) ? int(R__stl.size()) : 0; 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::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 &R__stl = trackPoints_; int R__n=(&R__stl) ? int(R__stl.size()) : 0; R__b << R__n; if(R__n) { std::vector::iterator R__k; for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) { R__b << (*R__k); } } } { std::map &R__stl = fitStatuses_; int R__n=(&R__stl) ? int(R__stl.size()) : 0; 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::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); } } } /* End of namespace genfit */