GBL/include/GblFitter.h

0001 /* Copyright 2013
0002  *   Authors: Sergey Yashchenko and Tadeas Bilka
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 /** @addtogroup genfit
0020  * @{
0021  */
0022
0023 #ifndef GblFitter_H
0024 #define GblFitter_H
0025
0026 #include "GblTrajectory.h"
0027 #include "AbsFitter.h"
0028 #include "AbsTrackRep.h"
0029 #include "GblFitterInfo.h"
0030 #include "GblFitStatus.h"
0031 #include "GblTrackSegmentController.h"
0032
0033 #include <map>
0034 #include <iostream>
0035
0036 #include <TMatrixD.h>
0037 #include <assert.h>
0038 #include <sstream>
0039
0040 #include <TMath.h>
0041 #include <TVector3.h>
0042
0043
0044 namespace genfit {
0045
0046   class GblTrackSegmentController;
0047
0048   /** @brief Generic GBL implementation
0049    *
0050    * The interface class to GBL track fit
0051    *
0052    */
0053   class GblFitter : public AbsFitter {
0054
0055   private:
0056     GblFitter(const GblFitter&);
0057     GblFitter& operator=(GblFitter const&);
0058
0059     std::string m_gblInternalIterations;
0060     bool m_enableScatterers;
0061     bool m_enableIntermediateScatterer;
0062     unsigned int m_externalIterations;
0063     unsigned int m_recalcJacobians;
0064
0065     // Minimum scattering sigma (will be squared and inverted...)
0066     double scatEpsilon;
0067     GblTrackSegmentController* m_segmentController;
0068
0069   public:
0070
0071     /**
0072      * Default (and only) constructor
0073      */
0074     GblFitter() : AbsFitter(), m_gblInternalIterations(""), m_enableScatterers(true), m_enableIntermediateScatterer(true), m_externalIterations(1), m_recalcJacobians(0), scatEpsilon(1.e-8), m_segmentController(nullptr) {;}
0075
0076     /**
0077      * Destructor
0078      */
0079     virtual ~GblFitter();
0080
0081     /**
0082      * @brief Set options of the fitter/GBL
0083      *
0084      * @return void
0085      * @param internalIterations GBL down-weighting in iterations. One letter (T,H,C) per iteration.
0086      *                           Seems "HH" is resonable for outliers. Default "" is OK.
0087      *                           Separate by comma for each external iteration
0088      *                           (usually last), e.g., for 3 iterations: ",,HH" to down-weight at last one
0089      * @param enableScatterers If false, no scatterers will be added to GBL trajectory
0090      * @param enableIntermediateScatterer True to simulate thick sctatterers by two thin scatterers
0091      *                                    1st at detector plane and intermediate between each two planes
0092      * @param externalIterations Sets number of times GblTrajectory.fit(...) will be called in processTrackWithRep(...).
0093      *                           Called external iterations. State is updated each time with GBL results.
0094      *                           If =0, GblFitterInfos will filled by reference states and GblFitStatus in track able
0095      *                           to construct simple GBL trajectory (for fit by GBL of output to Mille binary) / get list of points
0096      *                                    1st at detector plane and intermediate between each two planes
0097      * @param recalcJacobians Number of iteration up to which Jacobians should be recalculated / planes/meas updated after the fit.
0098      *                        0 = do not recalculate Jacobians. 1 = recalculate after first GBL fit. 2 = after 1st and 2nd GBL fit etc.
0099      */
0100     void setOptions(const std::string &internalIterations = "", bool enableScatterers = true, bool enableIntermediateScatterer = true, unsigned int externalIterations = 1, unsigned int recalcJacobians = 1) {
0101       m_externalIterations = externalIterations;
0102       m_gblInternalIterations = internalIterations;
0103       m_recalcJacobians = recalcJacobians;
0104       if (!enableScatterers)
0105         enableIntermediateScatterer = false;
0106       m_enableScatterers = enableScatterers;
0107       m_enableIntermediateScatterer = enableIntermediateScatterer;
0108     }
0109
0110
0111     /**
0112      * @brief Set multiple scattering options of the fitter/GBL
0113      *
0114      * @return void
0115      * @param enableScatterers If false, no scatterers will be added to GBL trajectory
0116      * @param enableIntermediateScatterer True to simulate thick sctatterers by two thin scatterers
0117      *                                    1st at detector plane and intermediate between each two planes
0118      */
0119     void setMSOptions(bool enableScatterers = true, bool enableIntermediateScatterer = true) {
0120       if (!enableScatterers)
0121         enableIntermediateScatterer = false;
0122       m_enableScatterers = enableScatterers;
0123       m_enableIntermediateScatterer = enableIntermediateScatterer;
0124     }
0125
0126     /**
0127      * @brief Evaluates moments of radiation length distribution from list of
0128      * material steps and computes parameters describing a corresponding thick scatterer.
0129      *
0130      * Based on input from Claus Kleinwort (DESY),
0131      * adapted for continuous material distribution represented by
0132      * a sum of step functions. Returned thick scatterer can be represented by two GBL scattering points.
0133      * Calculates variance of theta from total sum of radiation lengths
0134      * instead of summimg squares of individual deflection angle variances.
0135      *
0136      * @param length returned: Length of the track
0137      * @param theta returned: Variation of distribution of deflection angle
0138      * @param s returned: First moment of material scattering distribution
0139      * @param ds returned: Second moment (variance) of material scattering distribution
0140      * @param p Particle momentum magnitude (GeV/c)
0141      * @param mass Mass of particle (GeV/c/c)
0142      * @param steps Vector of material steps from (RKTrackRep) extrapolation
0143      * @return void
0144      */
0145     void getScattererFromMatList(double& length,
0146                                  double& theta, double& s, double& ds,
0147                                  const double p, const double mass, const double charge,
0148                                  const std::vector<genfit::MatStep>& steps) const;
0149
0150     /**
0151      * Performs fit on a Track.
0152      * Hit resorting currently supported (use only if necessary /wire chamber/ ... will
0153      * extrapolate along whole track to sort the hits).
0154      */
0155     void processTrackWithRep(Track* trk, const AbsTrackRep* rep, bool resortHits = false) override;
0156
0157     /**
0158      * @brief Propagate seed, populate track with scatterers
0159      * and GblFitterInfos with reference state set
0160      *
0161      * @param trk Track to attach with infos at given rep
0162      * @param rep TrackRep to which fitter info shall be attached
0163      * @return Length of track from extrapolations
0164      */
0165     double constructGblInfo(Track* trk, const AbsTrackRep* rep);
0166
0167     /**
0168      * @brief Populate all fitter infos in track for rep with
0169      * results of trajectory fit.
0170      *
0171      * Updates also seed state in track (from forward prediction
0172      * at first point)
0173      *
0174      * (The trajectory can only be cut before track end,
0175      * cannot have missing points in between (if valid))
0176      *
0177      *
0178      *
0179      * TODO ??
0180      * Re-construct all points in GblFitterInfos (updated) and collect them
0181      * in fit status
0182      *
0183      *
0184      * @param traj The fitted GblTrajectory
0185      * @param trk The track with fitter infos from whose points traj was created
0186      * @param rep The representation to which this fit status belong
0187      */
0188     void updateGblInfo(gbl::GblTrajectory& traj, genfit::Track* trk, const genfit::AbsTrackRep* rep);
0189
0190     /**
0191      * @brief Constructs all GBL points and returns them in vector
0192      * for trajectory construction
0193      *
0194      * @return const std::vector< gbl::GblPoint, std::allocator >&
0195      */
0196     std::vector<gbl::GblPoint> collectGblPoints(genfit::Track* trk, const genfit::AbsTrackRep* rep);
0197
0198     /**
0199      * @brief Remove all previous gbl fitter data from track
0200      * Also removes trackpoints without measurement
0201      *
0202      * @param trk The track
0203      * @param rep Representation for which to clean data
0204      * @return void
0205      */
0206     void cleanGblInfo(Track* trk, const AbsTrackRep* rep) const;
0207
0208     /**
0209      * @brief Sort hits in track by arc-len using extrapolation
0210      *
0211      * @param trk The track
0212      * @param rep The track representation
0213      * @return void
0214      */
0215     void sortHits(Track* trk, const AbsTrackRep* rep) const;
0216
0217     void setTrackSegmentController(GblTrackSegmentController* controler);
0218
0219
0220   public:
0221
0222     ClassDef(GblFitter, 2)
0223
0224   };
0225
0226 }  /* End of namespace genfit */
0227 /** @} */
0228
0229 #endif // GblFitter_H
0230