core/src/Track.cc

0001 /* Copyright 2008-2010, Technische Universitaet Muenchen,
0002    Authors: Christian Hoeppner & Sebastian Neubert & Johannes Rauch
0003
0004    This file is part of GENFIT.
0005
0006    GENFIT is free software: you can redistribute it and/or modify
0007    it under the terms of the GNU Lesser General Public License as published
0008    by the Free Software Foundation, either version 3 of the License, or
0009    (at your option) any later version.
0010
0011    GENFIT is distributed in the hope that it will be useful,
0012    but WITHOUT ANY WARRANTY; without even the implied warranty of
0013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
0014    GNU Lesser General Public License for more details.
0015
0016    You should have received a copy of the GNU Lesser General Public License
0017    along with GENFIT.  If not, see <http://www.gnu.org/licenses/>.
0018 */
0019
0020 #include "Track.h"
0021
0022 #include "Exception.h"
0023 #include "IO.h"
0024 #include "KalmanFitterInfo.h"
0025 #include "KalmanFitStatus.h"
0026 #include "PlanarMeasurement.h"
0027 #include "AbsMeasurement.h"
0028
0029 #include "WireTrackCandHit.h"
0030
0031 #include <algorithm>
0032 #include <map>
0033
0034 #include <TDatabasePDG.h>
0035 #include <TMath.h>
0036 #include <TBuffer.h>
0037
0038 //#include <glog/logging.h>
0039
0040 //#define DEBUG
0041
0042
0043 namespace genfit {
0044
0045 Track::Track() :
0046   cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(6), covSeed_(6)
0047 {
0048   ;
0049 }
0050
0051
0052 Track::Track(const TrackCand& trackCand, const MeasurementFactory<AbsMeasurement>& factory, AbsTrackRep* rep) :
0053   cardinalRep_(0), fitStatuses_(), stateSeed_(6), covSeed_(6)
0054 {
0055
0056   if (rep != nullptr)
0057     addTrackRep(rep);
0058
0059   createMeasurements(trackCand, factory);
0060
0061   // Copy seed information from candidate
0062   timeSeed_ = trackCand.getTimeSeed();
0063   stateSeed_ = trackCand.getStateSeed();
0064   covSeed_ = trackCand.getCovSeed();
0065
0066   mcTrackId_ = trackCand.getMcTrackId();
0067
0068   // fill cache
0069   fillPointsWithMeasurement();
0070
0071   checkConsistency();
0072 }
0073
0074 void
0075 Track::createMeasurements(const TrackCand& trackCand, const MeasurementFactory<AbsMeasurement>& factory)
0076 {
0077   // create the measurements using the factory.
0078   std::vector <AbsMeasurement*> factoryHits = factory.createMany(trackCand);
0079
0080   if (factoryHits.size() != trackCand.getNHits()) {
0081     Exception exc("Track::Track ==> factoryHits.size() != trackCand->getNHits()",__LINE__,__FILE__);
0082     exc.setFatal();
0083     throw exc;
0084   }
0085
0086   // create TrackPoints
0087   for (unsigned int i=0; i<factoryHits.size(); ++i){
0088     TrackPoint* tp = new TrackPoint(factoryHits[i], this);
0089     tp->setSortingParameter(trackCand.getHit(i)->getSortingParameter());
0090     insertPoint(tp);
0091   }
0092 }
0093
0094
0095 Track::Track(AbsTrackRep* trackRep, const TVectorD& stateSeed) :
0096   cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(stateSeed),
0097   covSeed_(TMatrixDSym::kUnit, TMatrixDSym(6))
0098 {
0099   addTrackRep(trackRep);
0100 }
0101
0102
0103 Track::Track(AbsTrackRep* trackRep, const TVector3& posSeed, const TVector3& momSeed) :
0104   cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(6),
0105   covSeed_(TMatrixDSym::kUnit, TMatrixDSym(6))
0106 {
0107   addTrackRep(trackRep);
0108   setStateSeed(posSeed, momSeed);
0109 }
0110
0111
0112 Track::Track(AbsTrackRep* trackRep, const TVectorD& stateSeed, const TMatrixDSym& covSeed) :
0113   cardinalRep_(0), fitStatuses_(), mcTrackId_(-1), timeSeed_(0), stateSeed_(stateSeed),
0114   covSeed_(covSeed)
0115 {
0116   addTrackRep(trackRep);
0117 }
0118
0119
0120 Track::Track(const Track& rhs) :
0121   TObject(rhs),
0122   cardinalRep_(rhs.cardinalRep_), mcTrackId_(rhs.mcTrackId_), timeSeed_(rhs.timeSeed_),
0123   stateSeed_(rhs.stateSeed_), covSeed_(rhs.covSeed_)
0124 {
0125   rhs.checkConsistency();
0126
0127   std::map<const AbsTrackRep*, AbsTrackRep*> oldRepNewRep;
0128
0129   for (std::vector<AbsTrackRep*>::const_iterator it=rhs.trackReps_.begin(); it!=rhs.trackReps_.end(); ++it) {
0130     AbsTrackRep* newRep = (*it)->clone();
0131     addTrackRep(newRep);
0132     oldRepNewRep[(*it)] = newRep;
0133   }
0134
0135   trackPoints_.reserve(rhs.trackPoints_.size());
0136   for (std::vector<TrackPoint*>::const_iterator it=rhs.trackPoints_.begin(); it!=rhs.trackPoints_.end(); ++it) {
0137     trackPoints_.push_back(new TrackPoint(**it, oldRepNewRep));
0138     trackPoints_.back()->setTrack(this);
0139   }
0140
0141   for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it=rhs.fitStatuses_.begin(); it!=rhs.fitStatuses_.end(); ++it) {
0142     setFitStatus(it->second->clone(), oldRepNewRep[it->first]);
0143   }
0144
0145   fillPointsWithMeasurement();
0146
0147   checkConsistency();
0148 }
0149
0150 Track& Track::operator=(Track other) {
0151   swap(other);
0152
0153   for (std::vector<TrackPoint*>::const_iterator it=trackPoints_.begin(); it!=trackPoints_.end(); ++it) {
0154     trackPoints_.back()->setTrack(this);
0155   }
0156
0157   fillPointsWithMeasurement();
0158
0159   checkConsistency();
0160
0161   return *this;
0162 }
0163
0164 void Track::swap(Track& other) {
0165   std::swap(this->trackReps_, other.trackReps_);
0166   std::swap(this->cardinalRep_, other.cardinalRep_);
0167   std::swap(this->trackPoints_, other.trackPoints_);
0168   std::swap(this->trackPointsWithMeasurement_, other.trackPointsWithMeasurement_);
0169   std::swap(this->fitStatuses_, other.fitStatuses_);
0170   std::swap(this->mcTrackId_, other.mcTrackId_);
0171   std::swap(this->timeSeed_, other.timeSeed_);
0172   std::swap(this->stateSeed_, other.stateSeed_);
0173   std::swap(this->covSeed_, other.covSeed_);
0174
0175 }
0176
0177 Track::~Track() {
0178   this->Clear();
0179 }
0180
0181 void Track::Clear(Option_t*)
0182 {
0183   // This function is needed for TClonesArray embedding.
0184   // FIXME: smarter containers or pointers needed ...
0185   for (size_t i = 0; i < trackPoints_.size(); ++i)
0186     delete trackPoints_[i];
0187
0188   trackPoints_.clear();
0189   trackPointsWithMeasurement_.clear();
0190
0191   for (std::map< const AbsTrackRep*, FitStatus* >::iterator it = fitStatuses_.begin(); it!= fitStatuses_.end(); ++it)
0192     delete it->second;
0193   fitStatuses_.clear();
0194
0195   for (size_t i = 0; i < trackReps_.size(); ++i)
0196     delete trackReps_[i];
0197   trackReps_.clear();
0198
0199   cardinalRep_ = 0;
0200
0201   mcTrackId_ = -1;
0202
0203   timeSeed_ = 0;
0204   stateSeed_.Zero();
0205   covSeed_.Zero();
0206 }
0207
0208
0209 TrackPoint* Track::getPoint(int id) const {
0210   if (id < 0)
0211     id += trackPoints_.size();
0212
0213   return trackPoints_.at(id);
0214 }
0215
0216
0217 TrackPoint* Track::getPointWithMeasurement(int id) const {
0218   if (id < 0)
0219     id += trackPointsWithMeasurement_.size();
0220
0221   return trackPointsWithMeasurement_.at(id);
0222 }
0223
0224
0225 TrackPoint* Track::getPointWithMeasurementAndFitterInfo(int id, const AbsTrackRep* rep) const {
0226   if (rep == nullptr)
0227     rep = getCardinalRep();
0228
0229   if (id >= 0) {
0230     int i = 0;
0231     for (std::vector<TrackPoint*>::const_iterator it = trackPointsWithMeasurement_.begin(); it != trackPointsWithMeasurement_.end(); ++it) {
0232       if ((*it)->hasFitterInfo(rep)) {
0233         if (id == i)
0234           return (*it);
0235         ++i;
0236       }
0237     }
0238   } else {
0239     // Search backwards.
0240     int i = -1;
0241     for (std::vector<TrackPoint*>::const_reverse_iterator it = trackPointsWithMeasurement_.rbegin(); it != trackPointsWithMeasurement_.rend(); ++it) {
0242       if ((*it)->hasFitterInfo(rep)) {
0243         if (id == i)
0244           return (*it);
0245         --i;
0246       }
0247     }
0248   }
0249
0250   // Not found, i.e. abs(id) > number of fitted TrackPoints
0251   return 0;
0252 }
0253
0254
0255 TrackPoint* Track::getPointWithFitterInfo(int id, const AbsTrackRep* rep) const {
0256   if (rep == nullptr)
0257     rep = getCardinalRep();
0258
0259   if (id >= 0) {
0260     int i = 0;
0261     for (std::vector<TrackPoint*>::const_iterator it = trackPoints_.begin(); it != trackPoints_.end(); ++it) {
0262       if ((*it)->hasFitterInfo(rep)) {
0263         if (id == i)
0264           return (*it);
0265         ++i;
0266       }
0267     }
0268   } else {
0269     // Search backwards.
0270     int i = -1;
0271     for (std::vector<TrackPoint*>::const_reverse_iterator it = trackPoints_.rbegin(); it != trackPoints_.rend(); ++it) {
0272       if ((*it)->hasFitterInfo(rep)) {
0273         if (id == i)
0274           return (*it);
0275         --i;
0276       }
0277     }
0278   }
0279
0280   // Not found, i.e. abs(id) > number of fitted TrackPoints
0281   return 0;
0282 }
0283
0284
0285 const MeasuredStateOnPlane& Track::getFittedState(int id, const AbsTrackRep* rep, bool biased) const {
0286   if (rep == nullptr)
0287     rep = getCardinalRep();
0288
0289   TrackPoint* point = getPointWithFitterInfo(id, rep);
0290   if (point == nullptr) {
0291     Exception exc("Track::getFittedState ==> no trackPoint with fitterInfo for rep",__LINE__,__FILE__);
0292     exc.setFatal();
0293     throw exc;
0294   }
0295   return point->getFitterInfo(rep)->getFittedState(biased);
0296 }
0297
0298
0299 int Track::getIdForRep(const AbsTrackRep* rep) const
0300 {
0301   for (size_t i = 0; i < trackReps_.size(); ++i)
0302     if (trackReps_[i] == rep)
0303       return i;
0304
0305   Exception exc("Track::getIdForRep ==> cannot find TrackRep in Track",__LINE__,__FILE__);
0306   exc.setFatal();
0307   throw exc;
0308 }
0309
0310
0311 bool Track::hasFitStatus(const AbsTrackRep* rep) const {
0312   if (rep == nullptr)
0313     rep = getCardinalRep();
0314
0315   if (fitStatuses_.find(rep) == fitStatuses_.end())
0316     return false;
0317
0318   return (fitStatuses_.at(rep) != nullptr);
0319 }
0320
0321
0322 bool Track::hasKalmanFitStatus(const AbsTrackRep* rep) const {
0323   if (rep == nullptr)
0324     rep = getCardinalRep();
0325
0326   if (fitStatuses_.find(rep) == fitStatuses_.end())
0327     return false;
0328
0329   return (dynamic_cast<KalmanFitStatus*>(fitStatuses_.at(rep)) != nullptr);
0330 }
0331
0332
0333 KalmanFitStatus* Track::getKalmanFitStatus(const AbsTrackRep* rep) const {
0334   return dynamic_cast<KalmanFitStatus*>(getFitStatus(rep));
0335 }
0336
0337
0338 void Track::setFitStatus(FitStatus* fitStatus, const AbsTrackRep* rep) {
0339   if (fitStatuses_.find(rep) != fitStatuses_.end())
0340     delete fitStatuses_.at(rep);
0341
0342   fitStatuses_[rep] = fitStatus;
0343 }
0344
0345
0346 void Track::setStateSeed(const TVector3& pos, const TVector3& mom) {
0347   stateSeed_.ResizeTo(6);
0348
0349   stateSeed_(0) = pos.X();
0350   stateSeed_(1) = pos.Y();
0351   stateSeed_(2) = pos.Z();
0352
0353   stateSeed_(3) = mom.X();
0354   stateSeed_(4) = mom.Y();
0355   stateSeed_(5) = mom.Z();
0356 }
0357
0358
0359
0360 void Track::insertPoint(TrackPoint* point, int id) {
0361
0362   point->setTrack(this);
0363
0364   #ifdef DEBUG
0365   debugOut << "Track::insertPoint at position " << id  << "\n";
0366   #endif
0367   assert(point!=nullptr);
0368   trackHasChanged();
0369
0370   point->setTrack(this);
0371
0372   if (trackPoints_.size() == 0) {
0373     trackPoints_.push_back(point);
0374
0375     if (point->hasRawMeasurements())
0376       trackPointsWithMeasurement_.push_back(point);
0377
0378     return;
0379   }
0380
0381   if (id == -1 || id == (int)trackPoints_.size()) {
0382     trackPoints_.push_back(point);
0383
0384     if (point->hasRawMeasurements())
0385       trackPointsWithMeasurement_.push_back(point);
0386
0387     deleteReferenceInfo(std::max(0, (int)trackPoints_.size()-2), (int)trackPoints_.size()-1);
0388
0389     // delete fitter infos if inserted point has a measurement
0390     if (point->hasRawMeasurements()) {
0391       deleteForwardInfo(-1, -1);
0392       deleteBackwardInfo(0, -2);
0393     }
0394
0395     return;
0396   }
0397
0398   // [-size, size-1] is the allowed range
0399   assert(id < (ssize_t)trackPoints_.size() || -id-1 <= (ssize_t)trackPoints_.size());
0400
0401   if (id < 0)
0402     id += trackPoints_.size() + 1;
0403
0404   // insert
0405   trackPoints_.insert(trackPoints_.begin() + id, point);  // insert inserts BEFORE
0406
0407   // delete fitter infos if inserted point has a measurement
0408   if (point->hasRawMeasurements()) {
0409     deleteForwardInfo(id, -1);
0410     deleteBackwardInfo(0, id);
0411   }
0412
0413   // delete reference info of neighbouring points
0414   deleteReferenceInfo(std::max(0, id-1), std::min((int)trackPoints_.size()-1, id+1));
0415
0416   fillPointsWithMeasurement();
0417 }
0418
0419
0420 void Track::insertPoints(std::vector<TrackPoint*> points, int id) {
0421
0422   int nBefore = getNumPoints();
0423   int n = points.size();
0424
0425   if (n == 0)
0426     return;
0427   if (n == 1) {
0428     insertPoint(points[0], id);
0429     return;
0430   }
0431
0432   for (std::vector<TrackPoint*>::iterator p = points.begin(); p != points.end(); ++p)
0433     (*p)->setTrack(this);
0434
0435   if (id == -1 || id == (int)trackPoints_.size()) {
0436     trackPoints_.insert(trackPoints_.end(), points.begin(), points.end());
0437
0438     deleteReferenceInfo(std::max(0, nBefore-1), nBefore);
0439
0440     deleteForwardInfo(nBefore, -1);
0441     deleteBackwardInfo(0, std::max(0, nBefore-1));
0442
0443     fillPointsWithMeasurement();
0444
0445     return;
0446   }
0447
0448
0449   assert(id < (ssize_t)trackPoints_.size() || -id-1 <= (ssize_t)trackPoints_.size());
0450
0451   if (id < 0)
0452     id += trackPoints_.size() + 1;
0453
0454
0455   // insert
0456   trackPoints_.insert(trackPoints_.begin() + id, points.begin(), points.end());  // insert inserts BEFORE
0457
0458   // delete fitter infos if inserted point has a measurement
0459   deleteForwardInfo(id, -1);
0460   deleteBackwardInfo(0, id+n);
0461
0462   // delete reference info of neighbouring points
0463   deleteReferenceInfo(std::max(0, id-1), std::min((int)trackPoints_.size()-1, id));
0464   deleteReferenceInfo(std::max(0, id+n-1), std::min((int)trackPoints_.size()-1, id+n));
0465
0466   fillPointsWithMeasurement();
0467 }
0468
0469
0470 void Track::deletePoint(int id) {
0471
0472   #ifdef DEBUG
0473   debugOut << "Track::deletePoint at position " << id  << "\n";
0474   #endif
0475
0476   trackHasChanged();
0477
0478   if (id < 0)
0479     id += trackPoints_.size();
0480   assert(id>0);
0481
0482
0483   // delete forwardInfo after point (backwardInfo before point) if deleted point has a measurement
0484   if (trackPoints_[id]->hasRawMeasurements()) {
0485     deleteForwardInfo(id, -1);
0486     deleteBackwardInfo(0, id-1);
0487   }
0488
0489   // delete reference info of neighbouring points
0490   deleteReferenceInfo(std::max(0, id-1), std::min((int)trackPoints_.size()-1, id+1));
0491
0492   // delete point
0493   std::vector<TrackPoint*>::iterator it = std::find(trackPointsWithMeasurement_.begin(), trackPointsWithMeasurement_.end(), trackPoints_[id]);
0494   if (it != trackPointsWithMeasurement_.end())
0495     trackPointsWithMeasurement_.erase(it);
0496
0497   delete trackPoints_[id];
0498   trackPoints_.erase(trackPoints_.begin()+id);
0499
0500   fillPointsWithMeasurement();
0501
0502 }
0503
0504
0505 void Track::insertMeasurement(AbsMeasurement* measurement, int id) {
0506   insertPoint(new TrackPoint(measurement, this), id);
0507 }
0508
0509 void Track::deleteFittedState(const genfit::AbsTrackRep* rep) {
0510   if(hasFitStatus(rep)) {
0511     delete fitStatuses_.at(rep);
0512     fitStatuses_.erase(rep);
0513   }
0514
0515   // delete FitterInfos related to the deleted TrackRep
0516   for (const auto& trackPoint : trackPoints_) {
0517     if(trackPoint->hasFitterInfo(rep)) {
0518       trackPoint->deleteFitterInfo(rep);
0519     }
0520   }
0521 }
0522
0523
0524 void Track::mergeTrack(const Track* other, int id) {
0525
0526   #ifdef DEBUG
0527   debugOut << "Track::mergeTrack\n";
0528   #endif
0529
0530   if (other->getNumPoints() == 0)
0531     return;
0532
0533   std::map<const AbsTrackRep*, AbsTrackRep*> otherRepThisRep;
0534   std::vector<const AbsTrackRep*> otherRepsToRemove;
0535
0536   for (std::vector<AbsTrackRep*>::const_iterator otherRep=other->trackReps_.begin(); otherRep!=other->trackReps_.end(); ++otherRep) {
0537     bool found(false);
0538     for (std::vector<AbsTrackRep*>::const_iterator thisRep=trackReps_.begin(); thisRep!=trackReps_.end(); ++thisRep) {
0539       if ((*thisRep)->isSame(*otherRep)) {
0540         otherRepThisRep[*otherRep] = *thisRep;
0541         #ifdef DEBUG
0542         debugOut << " map other rep " << *otherRep << " to " << (*thisRep) << "\n";
0543         #endif
0544         if (found) {
0545           Exception exc("Track::mergeTrack ==> more than one matching rep.",__LINE__,__FILE__);
0546           exc.setFatal();
0547           throw exc;
0548         }
0549         found = true;
0550         break;
0551       }
0552     }
0553     if (!found) {
0554       otherRepsToRemove.push_back(*otherRep);
0555       #ifdef DEBUG
0556       debugOut << " remove other rep " << *otherRep << "\n";
0557       #endif
0558     }
0559   }
0560
0561
0562   std::vector<TrackPoint*> points;
0563   points.reserve(other->getNumPoints());
0564
0565   for (std::vector<TrackPoint*>::const_iterator otherTp=other->trackPoints_.begin(); otherTp!=other->trackPoints_.end(); ++otherTp) {
0566     points.push_back(new TrackPoint(**otherTp, otherRepThisRep, &otherRepsToRemove));
0567   }
0568
0569   insertPoints(points, id);
0570 }
0571
0572
0573 void Track::addTrackRep(AbsTrackRep* trackRep) {
0574   trackReps_.push_back(trackRep);
0575   fitStatuses_[trackRep] = new FitStatus();
0576 }
0577
0578
0579 void Track::deleteTrackRep(int id) {
0580   if (id < 0)
0581     id += trackReps_.size();
0582
0583   AbsTrackRep* rep = trackReps_.at(id);
0584
0585   // update cardinalRep_
0586   if (int(cardinalRep_) == id)
0587     cardinalRep_ = 0; // reset
0588   else if (int(cardinalRep_) > id)
0589     --cardinalRep_; // make cardinalRep_ point to the same TrackRep before and after deletion
0590
0591   // delete FitterInfos related to the deleted TrackRep
0592   for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin(); pointIt != trackPoints_.end(); ++pointIt) {
0593     (*pointIt)->deleteFitterInfo(rep);
0594   }
0595
0596   // delete fitStatus
0597   delete fitStatuses_.at(rep);
0598   fitStatuses_.erase(rep);
0599
0600   // delete rep
0601   delete rep;
0602   trackReps_.erase(trackReps_.begin()+id);
0603 }
0604
0605
0606 void Track::setCardinalRep(int id) {
0607
0608   if (id < 0)
0609     id += trackReps_.size();
0610
0611   if (id >= 0 && (unsigned int)id < trackReps_.size())
0612     cardinalRep_ = id;
0613   else {
0614     cardinalRep_ = 0;
0615     errorOut << "Track::setCardinalRep: Attempted to set cardinalRep_ to a value out of bounds. Resetting  cardinalRep_ to 0." << std::endl;
0616   }
0617 }
0618
0619
0620 void Track::determineCardinalRep() {
0621
0622   // Todo: test
0623
0624   if (trackReps_.size() <= 1)
0625     return;
0626
0627   double minChi2(9.E99);
0628   const AbsTrackRep* bestRep(nullptr);
0629
0630   for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it = fitStatuses_.begin(); it != fitStatuses_.end(); ++it) {
0631     if (it->second->isFitConverged()) {
0632       if (it->second->getChi2() < minChi2) {
0633         minChi2 = it->second->getChi2();
0634         bestRep = it->first;
0635       }
0636     }
0637   }
0638
0639   if (bestRep != nullptr) {
0640     setCardinalRep(getIdForRep(bestRep));
0641   }
0642 }
0643
0644
0645 bool Track::sort() {
0646   #ifdef DEBUG
0647   debugOut << "Track::sort \n";
0648   #endif
0649
0650   int nPoints(trackPoints_.size());
0651   // original order
0652   std::vector<TrackPoint*> pointsBefore(trackPoints_);
0653
0654   // sort
0655   std::stable_sort(trackPoints_.begin(), trackPoints_.end(), TrackPointComparator());
0656
0657   // see where order changed
0658   int equalUntil(-1), equalFrom(nPoints);
0659   for (int i = 0; i<nPoints; ++i) {
0660     if (pointsBefore[i] == trackPoints_[i])
0661       equalUntil = i;
0662     else
0663       break;
0664   }
0665
0666   if (equalUntil == nPoints-1)
0667     return false; // sorting did not change anything
0668
0669
0670   trackHasChanged();
0671
0672   for (int i = nPoints-1; i>equalUntil; --i) {
0673     if (pointsBefore[i] == trackPoints_[i])
0674       equalFrom = i;
0675     else
0676       break;
0677   }
0678
0679   #ifdef DEBUG
0680   debugOut << "Track::sort. Equal up to (including) hit " << equalUntil << " and from (including) hit " << equalFrom << " \n";
0681   #endif
0682
0683   deleteForwardInfo(equalUntil+1, -1);
0684   deleteBackwardInfo(0, equalFrom-1);
0685
0686   deleteReferenceInfo(std::max(0, equalUntil+1), std::min((int)trackPoints_.size()-1, equalFrom-1));
0687
0688   fillPointsWithMeasurement();
0689
0690   return true;
0691 }
0692
0693
0694 bool Track::udpateSeed(int id, AbsTrackRep* rep, bool biased) {
0695   try {
0696     const MeasuredStateOnPlane& fittedState = getFittedState(id, rep, biased);
0697     setTimeSeed(fittedState.getTime());
0698     setStateSeed(fittedState.get6DState());
0699     setCovSeed(fittedState.get6DCov());
0700
0701     double fittedCharge = fittedState.getCharge();
0702
0703     for (unsigned int i = 0; i<trackReps_.size(); ++i) {
0704       if (trackReps_[i]->getPDGCharge() * fittedCharge < 0) {
0705         trackReps_[i]->switchPDGSign();
0706       }
0707     }
0708   }
0709   catch (Exception& e) {
0710     // in this case the original track seed will be used
0711     return false;
0712   }
0713   return true;
0714 }
0715
0716
0717 void Track::reverseTrackPoints() {
0718
0719   std::reverse(trackPoints_.begin(),trackPoints_.end());
0720
0721   deleteForwardInfo(0, -1);
0722   deleteBackwardInfo(0, -1);
0723   deleteReferenceInfo(0, -1);
0724
0725   fillPointsWithMeasurement();
0726 }
0727
0728
0729 void Track::switchPDGSigns(AbsTrackRep* rep) {
0730   if (rep != nullptr) {
0731     rep->switchPDGSign();
0732     return;
0733   }
0734
0735   for (unsigned int i = 0; i<trackReps_.size(); ++i) {
0736     trackReps_[i]->switchPDGSign();
0737   }
0738 }
0739
0740
0741 void Track::reverseTrack() {
0742   udpateSeed(-1); // set fitted state of last hit as new seed
0743   reverseMomSeed(); // flip momentum direction
0744   switchPDGSigns();
0745   reverseTrackPoints(); // also deletes all fitterInfos
0746 }
0747
0748
0749 void Track::deleteForwardInfo(int startId, int endId, const AbsTrackRep* rep) {
0750   #ifdef DEBUG
0751   debugOut << "Track::deleteForwardInfo from position " << startId  << " to " << endId << "\n";
0752   #endif
0753
0754   trackHasChanged();
0755
0756   if (startId < 0)
0757     startId += trackPoints_.size();
0758   if (endId < 0)
0759     endId += trackPoints_.size();
0760   endId += 1;
0761
0762   assert (endId >= startId);
0763
0764   for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) {
0765     if (rep != nullptr) {
0766       if ((*pointIt)->hasFitterInfo(rep))
0767         (*pointIt)->getFitterInfo(rep)->deleteForwardInfo();
0768     }
0769     else {
0770       const std::vector<AbsFitterInfo*> fitterInfos = (*pointIt)->getFitterInfos();
0771       for (std::vector<AbsFitterInfo*>::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end(); ++fitterInfoIt) {
0772         (*fitterInfoIt)->deleteForwardInfo();
0773       }
0774     }
0775   }
0776 }
0777
0778 void Track::deleteBackwardInfo(int startId, int endId, const AbsTrackRep* rep) {
0779
0780   #ifdef DEBUG
0781   debugOut << "Track::deleteBackwardInfo from position " << startId  << " to " << endId << "\n";
0782   #endif
0783
0784   trackHasChanged();
0785
0786   if (startId < 0)
0787     startId += trackPoints_.size();
0788   if (endId < 0)
0789     endId += trackPoints_.size();
0790   endId += 1;
0791
0792   assert (endId >= startId);
0793
0794
0795   for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) {
0796     if (rep != nullptr) {
0797       if ((*pointIt)->hasFitterInfo(rep))
0798         (*pointIt)->getFitterInfo(rep)->deleteBackwardInfo();
0799     }
0800     else {
0801       const std::vector<AbsFitterInfo*> fitterInfos = (*pointIt)->getFitterInfos();
0802       for (std::vector<AbsFitterInfo*>::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end(); ++fitterInfoIt) {
0803         (*fitterInfoIt)->deleteBackwardInfo();
0804       }
0805     }
0806   }
0807 }
0808
0809 void Track::deleteReferenceInfo(int startId, int endId, const AbsTrackRep* rep) {
0810
0811   #ifdef DEBUG
0812   debugOut << "Track::deleteReferenceInfo from position " << startId  << " to " << endId << "\n";
0813   #endif
0814
0815   trackHasChanged();
0816
0817   if (startId < 0)
0818     startId += trackPoints_.size();
0819   if (endId < 0)
0820     endId += trackPoints_.size();
0821   endId += 1;
0822
0823   assert (endId >= startId);
0824
0825   for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) {
0826     if (rep != nullptr) {
0827       if ((*pointIt)->hasFitterInfo(rep))
0828         (*pointIt)->getFitterInfo(rep)->deleteReferenceInfo();
0829     }
0830     else {
0831       std::vector<AbsFitterInfo*> fitterInfos = (*pointIt)->getFitterInfos();
0832       for (std::vector<AbsFitterInfo*>::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end(); ++fitterInfoIt) {
0833         (*fitterInfoIt)->deleteReferenceInfo();
0834       }
0835     }
0836   }
0837 }
0838
0839 void Track::deleteMeasurementInfo(int startId, int endId, const AbsTrackRep* rep) {
0840
0841   #ifdef DEBUG
0842   debugOut << "Track::deleteMeasurementInfo from position " << startId  << " to " << endId << "\n";
0843   #endif
0844
0845   trackHasChanged();
0846
0847   if (startId < 0)
0848     startId += trackPoints_.size();
0849   if (endId < 0)
0850     endId += trackPoints_.size();
0851   endId += 1;
0852
0853   assert (endId >= startId);
0854
0855   for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) {
0856     if (rep != nullptr) {
0857       if ((*pointIt)->hasFitterInfo(rep))
0858         (*pointIt)->getFitterInfo(rep)->deleteMeasurementInfo();
0859     }
0860     else {
0861       std::vector<AbsFitterInfo*> fitterInfos = (*pointIt)->getFitterInfos();
0862       for (std::vector<AbsFitterInfo*>::const_iterator fitterInfoIt = fitterInfos.begin(); fitterInfoIt != fitterInfos.end(); ++fitterInfoIt) {
0863         (*fitterInfoIt)->deleteMeasurementInfo();
0864       }
0865     }
0866   }
0867 }
0868
0869 void Track::deleteFitterInfo(int startId, int endId, const AbsTrackRep* rep) {
0870
0871   #ifdef DEBUG
0872   debugOut << "Track::deleteFitterInfo from position " << startId  << " to " << endId << "\n";
0873   #endif
0874
0875   trackHasChanged();
0876
0877   if (startId < 0)
0878     startId += trackPoints_.size();
0879   if (endId < 0)
0880     endId += trackPoints_.size();
0881   endId += 1;
0882
0883   assert (endId >= startId);
0884
0885   for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) {
0886     if (rep != nullptr) {
0887       if ((*pointIt)->hasFitterInfo(rep))
0888         (*pointIt)->deleteFitterInfo(rep);
0889     }
0890     else {
0891       for (std::vector<AbsTrackRep*>::const_iterator repIt = trackReps_.begin(); repIt != trackReps_.end(); ++repIt) {
0892         if ((*pointIt)->hasFitterInfo(*repIt))
0893           (*pointIt)->deleteFitterInfo(*repIt);
0894       }
0895     }
0896   }
0897 }
0898
0899
0900 double Track::getTrackLen(AbsTrackRep* rep, int startId, int endId) const {
0901
0902   if (startId < 0)
0903     startId += trackPoints_.size();
0904   if (endId < 0)
0905     endId += trackPoints_.size();
0906
0907   bool backwards(false);
0908   if (startId > endId) {
0909     double temp = startId;
0910     startId = endId;
0911     endId = temp;
0912     backwards = true;
0913   }
0914
0915   endId += 1;
0916
0917   if (rep == nullptr)
0918     rep = getCardinalRep();
0919
0920   double trackLen(0);
0921   StateOnPlane state;
0922
0923   for (std::vector<TrackPoint*>::const_iterator pointIt = trackPoints_.begin() + startId; pointIt != trackPoints_.begin() + endId; ++pointIt) {
0924     if (! (*pointIt)->hasFitterInfo(rep)) {
0925       Exception e("Track::getTracklength: trackPoint has no fitterInfo", __LINE__,__FILE__);
0926       throw e;
0927     }
0928
0929     if (pointIt != trackPoints_.begin() + startId) {
0930       trackLen += rep->extrapolateToPlane(state, (*pointIt)->getFitterInfo(rep)->getPlane());
0931     }
0932
0933     state = (*pointIt)->getFitterInfo(rep)->getFittedState();
0934   }
0935
0936   if (backwards)
0937     trackLen *= -1.;
0938
0939   return trackLen;
0940 }
0941
0942
0943 TrackCand* Track::constructTrackCand() const {
0944   TrackCand* cand = new TrackCand();
0945
0946   cand->setTime6DSeedAndPdgCode(timeSeed_, stateSeed_, getCardinalRep()->getPDG());
0947   cand->setCovSeed(covSeed_);
0948   cand->setMcTrackId(mcTrackId_);
0949
0950   for (unsigned int i = 0; i < trackPointsWithMeasurement_.size(); ++i) {
0951     const TrackPoint* tp = trackPointsWithMeasurement_[i];
0952     const std::vector< AbsMeasurement* >& measurements = tp->getRawMeasurements();
0953
0954     for (unsigned int j = 0; j < measurements.size(); ++j) {
0955       const AbsMeasurement* m = measurements[j];
0956       TrackCandHit* tch;
0957
0958       int planeId = -1;
0959       if (dynamic_cast<const PlanarMeasurement*>(m)) {
0960         planeId = static_cast<const PlanarMeasurement*>(m)->getPlaneId();
0961       }
0962
0963       if (m->isLeftRightMeasurement()) {
0964         tch = new WireTrackCandHit(m->getDetId(), m->getHitId(), planeId,
0965             tp->getSortingParameter(), m->getLeftRightResolution());
0966       }
0967       else {
0968         tch = new TrackCandHit(m->getDetId(), m->getHitId(), planeId,
0969             tp->getSortingParameter());
0970       }
0971       cand->addHit(tch);
0972     }
0973   }
0974
0975   return cand;
0976 }
0977
0978
0979 double Track::getTOF(AbsTrackRep* rep, int startId, int endId) const {
0980
0981   if (startId < 0)
0982     startId += trackPoints_.size();
0983   if (endId < 0)
0984     endId += trackPoints_.size();
0985
0986   if (startId > endId) {
0987     std::swap(startId, endId);
0988   }
0989
0990   endId += 1;
0991
0992   if (rep == nullptr)
0993     rep = getCardinalRep();
0994
0995   StateOnPlane state;
0996
0997   const TrackPoint* startPoint(trackPoints_[startId]);
0998   const TrackPoint* endPoint(trackPoints_[endId]);
0999
1000   if (!startPoint->hasFitterInfo(rep)
1001       || !endPoint->hasFitterInfo(rep)) {
1002       Exception e("Track::getTOF: trackPoint has no fitterInfo", __LINE__,__FILE__);
1003       throw e;
1004     }
1005
1006   double tof = (rep->getTime(endPoint->getFitterInfo(rep)->getFittedState())
1007         - rep->getTime(startPoint->getFitterInfo(rep)->getFittedState()));
1008   return tof;
1009 }
1010
1011
1012 void Track::fixWeights(AbsTrackRep* rep, int startId, int endId) {
1013
1014   if (startId < 0)
1015     startId += trackPoints_.size();
1016   if (endId < 0)
1017     endId += trackPoints_.size();
1018
1019   assert(startId >= 0);
1020   assert(startId <= endId);
1021   assert(endId <= (int)trackPoints_.size());
1022
1023   std::vector< AbsFitterInfo* > fis;
1024
1025   for (std::vector<TrackPoint*>::iterator tp = trackPoints_.begin() + startId; tp != trackPoints_.begin() + endId; ++tp) {
1026     fis.clear();
1027     if (rep == nullptr) {
1028       fis = (*tp)->getFitterInfos();
1029     }
1030     else if ((*tp)->hasFitterInfo(rep)) {
1031       fis.push_back((*tp)->getFitterInfo(rep));
1032     }
1033
1034     for (std::vector< AbsFitterInfo* >::iterator fi = fis.begin(); fi != fis.end(); ++fi) {
1035       KalmanFitterInfo* kfi = dynamic_cast<KalmanFitterInfo*>(*fi);
1036       if (kfi == nullptr)
1037         continue;
1038
1039       kfi->fixWeights();
1040     }
1041   }
1042 }
1043
1044
1045 void Track::prune(const Option_t* option) {
1046
1047   PruneFlags f;
1048   f.setFlags(option);
1049
1050   for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it=fitStatuses_.begin(); it!=fitStatuses_.end(); ++it) {
1051     it->second->getPruneFlags().setFlags(option);
1052   }
1053
1054   // prune trackPoints
1055   if (f.hasFlags("F") || f.hasFlags("L")) {
1056     TrackPoint* firstPoint = getPointWithFitterInfo(0);
1057     TrackPoint* lastPoint = getPointWithFitterInfo(-1);
1058     for (unsigned int i = 0; i<trackPoints_.size(); ++i) {
1059       if (trackPoints_[i] == firstPoint && f.hasFlags("F"))
1060         continue;
1061
1062       if (trackPoints_[i] == lastPoint && f.hasFlags("L"))
1063         continue;
1064
1065       delete trackPoints_[i];
1066       trackPoints_.erase(trackPoints_.begin()+i);
1067       --i;
1068     }
1069   }
1070
1071   // prune TrackReps
1072   if (f.hasFlags("C")) {
1073     for (unsigned int i = 0; i < trackReps_.size(); ++i) {
1074       if (i != cardinalRep_) {
1075         deleteTrackRep(i);
1076         --i;
1077       }
1078     }
1079   }
1080
1081
1082   // from remaining trackPoints: prune measurementsOnPlane, unneeded fitterInfoStuff
1083   for (unsigned int i = 0; i<trackPoints_.size(); ++i) {
1084     if (f.hasFlags("W"))
1085       trackPoints_[i]->deleteRawMeasurements();
1086
1087     std::vector< AbsFitterInfo* > fis =  trackPoints_[i]->getFitterInfos();
1088     for (unsigned int j = 0; j<fis.size(); ++j) {
1089
1090       if (i == 0 && f.hasFlags("FLI"))
1091         fis[j]->deleteForwardInfo();
1092       else if (i == trackPoints_.size()-1 && f.hasFlags("FLI"))
1093         fis[j]->deleteBackwardInfo();
1094       else if (f.hasFlags("FI"))
1095         fis[j]->deleteForwardInfo();
1096       else if (f.hasFlags("LI"))
1097         fis[j]->deleteBackwardInfo();
1098
1099       if (f.hasFlags("U") && dynamic_cast<KalmanFitterInfo*>(fis[j]) != nullptr) {
1100         static_cast<KalmanFitterInfo*>(fis[j])->deletePredictions();
1101       }
1102
1103       // also delete reference info if points have been removed since it is invalid then!
1104       if (f.hasFlags("R") or f.hasFlags("F") or f.hasFlags("L"))
1105         fis[j]->deleteReferenceInfo();
1106       if (f.hasFlags("M"))
1107         fis[j]->deleteMeasurementInfo();
1108     }
1109   }
1110
1111   fillPointsWithMeasurement();
1112
1113   #ifdef DEBUG
1114   debugOut << "pruned Track: "; Print();
1115   #endif
1116
1117 }
1118
1119
1120 void Track::Print(const Option_t* option) const {
1121   TString opt = option;
1122   opt.ToUpper();
1123   if (opt.Contains("C")) { // compact
1124
1125     printOut << "\n    ";
1126     for (unsigned int i=0; i<trackPoints_.size(); ++i) {
1127
1128       int color = 32*(size_t)(trackPoints_[i]) % 15;
1129       switch (color) {
1130         case 0:
1131           printOut<<"\033[1;30m";
1132           break;
1133         case 1:
1134           printOut<<"\033[0;34m";
1135           break;
1136         case 2:
1137           printOut<<"\033[1;34m";
1138           break;
1139         case 3:
1140           printOut<<"\033[0;32m";
1141           break;
1142         case 4:
1143           printOut<<"\033[1;32m";
1144           break;
1145         case 5:
1146           printOut<<"\033[0;36m";
1147           break;
1148         case 6:
1149           printOut<<"\033[1;36m";
1150           break;
1151         case 7:
1152           printOut<<"\033[0;31m";
1153           break;
1154         case 8:
1155           printOut<<"\033[1;31m";
1156           break;
1157         case 9:
1158           printOut<<"\033[0;35m";
1159           break;
1160         case 10:
1161           printOut<<"\033[1;35m";
1162           break;
1163         case 11:
1164           printOut<<"\033[0;33m";
1165           break;
1166         case 12:
1167           printOut<<"\033[1;33m";
1168           break;
1169         case 13:
1170           printOut<<"\033[0;37m";
1171           break;
1172         default:
1173           ;
1174       }
1175       printOut << trackPoints_[i] << "\033[00m  ";
1176     }
1177     printOut << "\n";
1178
1179     printOut << "   ";
1180     for (unsigned int i=0; i<trackPoints_.size(); ++i) {
1181       printf("% -9.3g  ", trackPoints_[i]->getSortingParameter());
1182     }
1183
1184     for (std::vector<AbsTrackRep*>::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) {
1185       printOut << "\n" << getIdForRep(*rep) << "   ";
1186       for (unsigned int i=0; i<trackPoints_.size(); ++i) {
1187         if (! trackPoints_[i]->hasFitterInfo(*rep)) {
1188           printOut << "           ";
1189           continue;
1190         }
1191         AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep);
1192         if (fi->hasMeasurements())
1193           printOut << "M";
1194         else
1195           printOut << " ";
1196
1197         if (fi->hasReferenceState())
1198           printOut << "R";
1199         else
1200           printOut << " ";
1201
1202         printOut << "         ";
1203       }
1204       printOut << "\n";
1205
1206       printOut << " -> ";
1207       for (unsigned int i=0; i<trackPoints_.size(); ++i) {
1208         if (! trackPoints_[i]->hasFitterInfo(*rep)) {
1209           printOut << "           ";
1210           continue;
1211         }
1212         AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep);
1213         if (fi->hasForwardPrediction())
1214           printOut << "P";
1215         else
1216           printOut << " ";
1217
1218         if (fi->hasForwardUpdate())
1219           printOut << "U";
1220         else
1221           printOut << " ";
1222
1223         printOut << "         ";
1224       }
1225       printOut << "\n";
1226
1227       printOut << " <- ";
1228       for (unsigned int i=0; i<trackPoints_.size(); ++i) {
1229         if (! trackPoints_[i]->hasFitterInfo(*rep)) {
1230           printOut << "           ";
1231           continue;
1232         }
1233         AbsFitterInfo* fi = trackPoints_[i]->getFitterInfo(*rep);
1234         if (fi->hasBackwardPrediction())
1235           printOut << "P";
1236         else
1237           printOut << " ";
1238
1239         if (fi->hasBackwardUpdate())
1240           printOut << "U";
1241         else
1242           printOut << " ";
1243
1244         printOut << "         ";
1245       }
1246
1247       printOut << "\n";
1248
1249     } //end loop over reps
1250
1251     printOut << "\n";
1252     return;
1253   }
1254
1255
1256
1257   printOut << "=======================================================================================\n";
1258   printOut << "genfit::Track, containing " << trackPoints_.size() << " TrackPoints and " << trackReps_.size() << " TrackReps.\n";
1259   printOut << " Seed state: "; stateSeed_.Print();
1260
1261   for (unsigned int i=0; i<trackReps_.size(); ++i) {
1262     printOut << " TrackRep Nr. " << i;
1263     if (i == cardinalRep_)
1264       printOut << " (This is the cardinal rep)";
1265     printOut << "\n";
1266     trackReps_[i]->Print();
1267   }
1268
1269   printOut << "---------------------------------------------------------------------------------------\n";
1270
1271   for (unsigned int i=0; i<trackPoints_.size(); ++i) {
1272     printOut << "TrackPoint Nr. " << i << "\n";
1273     trackPoints_[i]->Print();
1274     printOut << "..........................................................................\n";
1275   }
1276
1277   for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it=fitStatuses_.begin(); it!=fitStatuses_.end(); ++it) {
1278     it->second->Print();
1279   }
1280
1281   printOut << "=======================================================================================\n";
1282
1283 }
1284
1285
1286 void Track::checkConsistency() const {
1287
1288   // cppcheck-suppress unreadVariable
1289   bool consistent = true;
1290   std::stringstream failures;
1291
1292   std::map<const AbsTrackRep*, const KalmanFitterInfo*> prevFis;
1293
1294   // check if seed is 6D
1295   if (stateSeed_.GetNrows() != 6) {
1296     failures << "Track::checkConsistency(): stateSeed_ dimension != 6" << std::endl;
1297     // cppcheck-suppress unreadVariable
1298     consistent = false;
1299   }
1300
1301   if (covSeed_.GetNrows() != 6) {
1302     failures << "Track::checkConsistency(): covSeed_ dimension != 6" << std::endl;
1303     // cppcheck-suppress unreadVariable
1304     consistent = false;
1305   }
1306
1307   if (covSeed_.Max() == 0.) {
1308     // Nota bene: The consistency is not set to false when this occurs, because it does not break the consistency of
1309     // the track. However, when something else fails we keep this as additional error information.
1310     failures << "Track::checkConsistency(): Warning: covSeed_ is zero" << std::endl;
1311   }
1312
1313   // check if correct number of fitStatuses
1314   if (fitStatuses_.size() != trackReps_.size()) {
1315     failures << "Track::checkConsistency(): Number of fitStatuses is != number of TrackReps " << std::endl;
1316     // cppcheck-suppress unreadVariable
1317     consistent = false;
1318   }
1319
1320   // check if cardinalRep_ is in range of trackReps_
1321   if (trackReps_.size() && cardinalRep_ >= trackReps_.size()) {
1322     failures << "Track::checkConsistency(): cardinalRep id " << cardinalRep_ << " out of bounds" << std::endl;
1323     // cppcheck-suppress unreadVariable
1324     consistent = false;
1325   }
1326
1327   for (std::vector<AbsTrackRep*>::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) {
1328     // check for nullptr
1329     if ((*rep) == nullptr) {
1330       failures << "Track::checkConsistency(): TrackRep is nullptr" << std::endl;
1331       // cppcheck-suppress unreadVariable
1332       consistent = false;
1333     }
1334
1335     // check for valid pdg code
1336     TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle((*rep)->getPDG());
1337     if (particle == nullptr) {
1338       failures << "Track::checkConsistency(): TrackRep pdg ID " << (*rep)->getPDG() << " is not valid" << std::endl;
1339       // cppcheck-suppress unreadVariable
1340       consistent = false;
1341     }
1342
1343     // check if corresponding FitStatus is there
1344     if (fitStatuses_.find(*rep) == fitStatuses_.end() and fitStatuses_.find(*rep)->second != nullptr) {
1345       failures << "Track::checkConsistency(): No FitStatus for Rep or FitStatus is nullptr" << std::endl;
1346       // cppcheck-suppress unreadVariable
1347       consistent = false;
1348     }
1349   }
1350
1351   // check TrackPoints
1352   for (std::vector<TrackPoint*>::const_iterator tp = trackPoints_.begin(); tp != trackPoints_.end(); ++tp) {
1353     // check for nullptr
1354     if ((*tp) == nullptr) {
1355       failures << "Track::checkConsistency(): TrackPoint is nullptr" << std::endl;
1356       // cppcheck-suppress unreadVariable
1357       consistent = false;
1358     }
1359     // check if trackPoint points back to this track
1360     if ((*tp)->getTrack() != this) {
1361       failures << "Track::checkConsistency(): TrackPoint does not point back to this track" << std::endl;
1362       // cppcheck-suppress unreadVariable
1363       consistent = false;
1364     }
1365
1366     // check rawMeasurements
1367     const std::vector<AbsMeasurement*>& rawMeasurements = (*tp)->getRawMeasurements();
1368     for (std::vector<AbsMeasurement*>::const_iterator m = rawMeasurements.begin(); m != rawMeasurements.end(); ++m) {
1369       // check for nullptr
1370       if ((*m) == nullptr) {
1371         failures << "Track::checkConsistency(): Measurement is nullptr" << std::endl;
1372     // cppcheck-suppress unreadVariable
1373         consistent = false;
1374       }
1375       // check if measurement points back to TrackPoint
1376       if ((*m)->getTrackPoint() != *tp) {
1377         failures << "Track::checkConsistency(): Measurement does not point back to correct TrackPoint" << std::endl;
1378     // cppcheck-suppress unreadVariable
1379         consistent = false;
1380       }
1381     }
1382
1383     // check fitterInfos
1384     std::vector<AbsFitterInfo*> fitterInfos = (*tp)->getFitterInfos();
1385     for (std::vector<AbsFitterInfo*>::const_iterator fi = fitterInfos.begin(); fi != fitterInfos.end(); ++fi) {
1386       // check for nullptr
1387       if ((*fi) == nullptr) {
1388         failures << "Track::checkConsistency(): FitterInfo is nullptr. TrackPoint: " << *tp << std::endl;
1389     // cppcheck-suppress unreadVariable
1390         consistent = false;
1391       }
1392
1393       // check if fitterInfos point to valid TrackReps in trackReps_
1394       int mycount (0);
1395       for (std::vector<AbsTrackRep*>::const_iterator rep = trackReps_.begin(); rep != trackReps_.end(); ++rep) {
1396         if ( (*rep) == (*fi)->getRep() ) {
1397           ++mycount;
1398         }
1399       }
1400       if (mycount ==  0) {
1401         failures << "Track::checkConsistency(): fitterInfo points to TrackRep which is not in Track" << std::endl;
1402     // cppcheck-suppress unreadVariable
1403         consistent = false;
1404       }
1405
1406       if (!( (*fi)->checkConsistency(&(this->getFitStatus((*fi)->getRep())->getPruneFlags())) ) ) {
1407         failures << "Track::checkConsistency(): FitterInfo not consistent. TrackPoint: " << *tp << std::endl;
1408     // cppcheck-suppress unreadVariable
1409         consistent = false;
1410       }
1411
1412       if (dynamic_cast<KalmanFitterInfo*>(*fi) != nullptr) {
1413         if (prevFis[(*fi)->getRep()] != nullptr &&
1414             static_cast<KalmanFitterInfo*>(*fi)->hasReferenceState() &&
1415             prevFis[(*fi)->getRep()]->hasReferenceState() ) {
1416           double len = static_cast<KalmanFitterInfo*>(*fi)->getReferenceState()->getForwardSegmentLength();
1417           double prevLen = prevFis[(*fi)->getRep()]->getReferenceState()->getBackwardSegmentLength();
1418           if (fabs(prevLen + len) > 1E-10 ) {
1419             failures << "Track::checkConsistency(): segment lengths of reference states for rep " << (*fi)->getRep() << " (id " << getIdForRep((*fi)->getRep()) << ") at TrackPoint " << (*tp) << " don't match" << std::endl;
1420             failures << prevLen << " + " << len << " = " << prevLen + len << std::endl;
1421             failures << "TrackPoint " << *tp << ", FitterInfo " << *fi << ", rep " << getIdForRep((*fi)->getRep()) << std::endl;
1422         // cppcheck-suppress unreadVariable
1423             consistent = false;
1424           }
1425         }
1426
1427         prevFis[(*fi)->getRep()] = static_cast<KalmanFitterInfo*>(*fi);
1428       }
1429       else
1430         prevFis[(*fi)->getRep()] = nullptr;
1431
1432     } // end loop over FitterInfos
1433
1434   } // end loop over TrackPoints
1435
1436
1437   // check trackPointsWithMeasurement_
1438   std::vector<TrackPoint*> trackPointsWithMeasurement;
1439
1440   for (std::vector<TrackPoint*>::const_iterator it = trackPoints_.begin(); it != trackPoints_.end(); ++it) {
1441     if ((*it)->hasRawMeasurements()) {
1442       trackPointsWithMeasurement.push_back(*it);
1443     }
1444   }
1445
1446   if (trackPointsWithMeasurement.size() != trackPointsWithMeasurement_.size()) {
1447     failures << "Track::checkConsistency(): trackPointsWithMeasurement_ has incorrect size" << std::endl;
1448     // cppcheck-suppress unreadVariable
1449     consistent = false;
1450   }
1451
1452   for (unsigned int i = 0; i < trackPointsWithMeasurement.size(); ++i) {
1453     if (trackPointsWithMeasurement[i] != trackPointsWithMeasurement_[i]) {
1454       failures << "Track::checkConsistency(): trackPointsWithMeasurement_ is not correct" << std::endl;
1455       failures << "has         id " << i << ", address " << trackPointsWithMeasurement_[i] << std::endl;
1456       failures << "should have id " << i << ", address " << trackPointsWithMeasurement[i] << std::endl;
1457       // cppcheck-suppress unreadVariable
1458       consistent = false;
1459     }
1460   }
1461
1462   if (not consistent) {
1463     throw genfit::Exception(failures.str(), __LINE__, __FILE__);
1464   }
1465 }
1466
1467
1468 void Track::trackHasChanged() {
1469
1470   #ifdef DEBUG
1471   debugOut << "Track::trackHasChanged \n";
1472   #endif
1473
1474   if (fitStatuses_.empty())
1475     return;
1476
1477   for (std::map< const AbsTrackRep*, FitStatus* >::const_iterator it=fitStatuses_.begin(); it!=fitStatuses_.end(); ++it) {
1478     it->second->setHasTrackChanged();
1479   }
1480 }
1481
1482
1483 void Track::fillPointsWithMeasurement() {
1484   trackPointsWithMeasurement_.clear();
1485   trackPointsWithMeasurement_.reserve(trackPoints_.size());
1486
1487   for (std::vector<TrackPoint*>::const_iterator it = trackPoints_.begin(); it != trackPoints_.end(); ++it) {
1488     if ((*it)->hasRawMeasurements()) {
1489       trackPointsWithMeasurement_.push_back(*it);
1490     }
1491   }
1492 }
1493
1494
1495 void Track::Streamer(TBuffer &R__b)
1496 {
1497    // Stream an object of class genfit::Track.
1498   const bool streamTrackPoints = true; // debugging option
1499    //This works around a msvc bug and should be harmless on other platforms
1500    typedef ::genfit::Track thisClass;
1501    UInt_t R__s, R__c;
1502    if (R__b.IsReading()) {
1503       this->Clear();
1504       Version_t R__v = R__b.ReadVersion(&R__s, &R__c); if (R__v) { }
1505       TObject::Streamer(R__b);
1506       {
1507         std::vector<AbsTrackRep*> &R__stl =  trackReps_;
1508         R__stl.clear();
1509         TClass *R__tcl1 = TBuffer::GetClass(typeid(genfit::AbsTrackRep));
1510         if (R__tcl1==0) {
1511           Error("trackReps_ streamer","Missing the TClass object for genfit::AbsTrackRep!");
1512           return;
1513         }
1514         int R__i, R__n;
1515         R__b >> R__n;
1516         R__stl.reserve(R__n);
1517         for (R__i = 0; R__i < R__n; R__i++) {
1518           genfit::AbsTrackRep* R__t;
1519           R__b >> R__t;
1520           R__stl.push_back(R__t);
1521         }
1522       }
1523       R__b >> cardinalRep_;
1524       if (streamTrackPoints)
1525       {
1526         std::vector<TrackPoint*> &R__stl =  trackPoints_;
1527         R__stl.clear();
1528         TClass *R__tcl1 = TBuffer::GetClass(typeid(genfit::TrackPoint));
1529         if (R__tcl1==0) {
1530           Error("trackPoints_ streamer","Missing the TClass object for genfit::TrackPoint!");
1531           return;
1532         }
1533         int R__i, R__n;
1534         R__b >> R__n;
1535         R__stl.reserve(R__n);
1536         for (R__i = 0; R__i < R__n; R__i++) {
1537           genfit::TrackPoint* R__t;
1538           R__t = (genfit::TrackPoint*)R__b.ReadObjectAny(R__tcl1);
1539           R__t->setTrack(this);
1540           R__t->fixupRepsForReading();
1541           R__stl.push_back(R__t);
1542         }
1543       }
1544       {
1545         std::map<const AbsTrackRep*,FitStatus*> &R__stl =  fitStatuses_;
1546         R__stl.clear();
1547         TClass *R__tcl1 = TBuffer::GetClass(typeid(genfit::AbsTrackRep));
1548         if (R__tcl1==0) {
1549           Error("fitStatuses_ streamer","Missing the TClass object for genfit::AbsTrackRep!");
1550           return;
1551         }
1552         TClass *R__tcl2 = TBuffer::GetClass(typeid(genfit::FitStatus));
1553         if (R__tcl2==0) {
1554           Error("fitStatuses_ streamer","Missing the TClass object for genfit::FitStatus!");
1555           return;
1556         }
1557         int R__i, R__n;
1558         R__b >> R__n;
1559         for (R__i = 0; R__i < R__n; R__i++) {
1560           int id;
1561           R__b >> id;
1562           genfit::FitStatus* R__t2;
1563           R__t2 = (genfit::FitStatus*)R__b.ReadObjectAny(R__tcl2);
1564
1565           R__stl[this->getTrackRep(id)] = R__t2;
1566         }
1567       }
1568       // timeSeed_ was introduced in version 3.  If reading an earlier
1569       // version, default to zero.
1570       if (R__v >= 3) {
1571     R__b >> timeSeed_;
1572       } else {
1573     timeSeed_ = 0;
1574       }
1575       stateSeed_.Streamer(R__b);
1576       covSeed_.Streamer(R__b);
1577       R__b.CheckByteCount(R__s, R__c, thisClass::IsA());
1578
1579       fillPointsWithMeasurement();
1580    } else {
1581       R__c = R__b.WriteVersion(thisClass::IsA(), kTRUE);
1582       TObject::Streamer(R__b);
1583       {
1584         std::vector<AbsTrackRep*> &R__stl =  trackReps_;
1585         int R__n=int(R__stl.size());
1586         R__b << R__n;
1587         if(R__n) {
1588           TClass *R__tcl1 = TBuffer::GetClass(typeid(genfit::AbsTrackRep));
1589           if (R__tcl1==0) {
1590             Error("trackReps_ streamer","Missing the TClass object for genfit::AbsTrackRep!");
1591             return;
1592           }
1593           std::vector<AbsTrackRep*>::iterator R__k;
1594           for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1595             R__b << *R__k;
1596           }
1597         }
1598       }
1599       R__b << cardinalRep_;
1600       if (streamTrackPoints)
1601       {
1602         std::vector<TrackPoint*> &R__stl =  trackPoints_;
1603         int R__n=int(R__stl.size());
1604         R__b << R__n;
1605         if(R__n) {
1606           std::vector<TrackPoint*>::iterator R__k;
1607           for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1608             R__b << (*R__k);
1609           }
1610         }
1611       }
1612       {
1613         std::map<const AbsTrackRep*,FitStatus*> &R__stl =  fitStatuses_;
1614         int R__n=int(R__stl.size());
1615         R__b << R__n;
1616         if(R__n) {
1617           TClass *R__tcl1 = TBuffer::GetClass(typeid(genfit::AbsTrackRep));
1618           if (R__tcl1==0) {
1619             Error("fitStatuses_ streamer","Missing the TClass object for genfit::AbsTrackRep!");
1620             return;
1621           }
1622           TClass *R__tcl2 = TBuffer::GetClass(typeid(genfit::FitStatus));
1623           if (R__tcl2==0) {
1624             Error("fitStatuses_ streamer","Missing the TClass object for genfit::FitStatus!");
1625             return;
1626           }
1627           std::map<const AbsTrackRep*,FitStatus*>::iterator R__k;
1628           for (R__k = R__stl.begin(); R__k != R__stl.end(); ++R__k) {
1629             int id = this->getIdForRep((*R__k).first);
1630             R__b << id;
1631             R__b << ((*R__k).second);
1632           }
1633         }
1634       }
1635       R__b << timeSeed_;
1636       stateSeed_.Streamer(R__b);
1637       covSeed_.Streamer(R__b);
1638       R__b.SetByteCount(R__c, kTRUE);
1639    }
1640 }
1641
1642 void Track::deleteTrackPointsAndFitStatus() {
1643   for (size_t i = 0; i < trackPoints_.size(); ++i)
1644     delete trackPoints_[i];
1645
1646   trackPoints_.clear();
1647   trackPointsWithMeasurement_.clear();
1648
1649   for (std::map< const AbsTrackRep*, FitStatus* >::iterator it = fitStatuses_.begin(); it!= fitStatuses_.end(); ++it)
1650     delete it->second;
1651   fitStatuses_.clear();
1652 }
1653
1654 } /* End of namespace genfit */