Commit 10a2116a authored by Simon Spannagel's avatar Simon Spannagel
Browse files

Add first version of Millepede algorithm, taken from LHCb's Kepler framework

parent 56c6a269
Pipeline #321837 passed with stages
in 4 minutes and 20 seconds
# Define module and return the generated name as ALGORITHM_NAME
CORRYVRECKAN_ALGORITHM(ALGORITHM_NAME)
# Add source files to library
CORRYVRECKAN_ALGORITHM_SOURCES(${ALGORITHM_NAME}
Millepede.cpp
# ADD SOURCE FILES HERE...
)
# Provide standard install target
CORRYVRECKAN_ALGORITHM_INSTALL(${ALGORITHM_NAME})
This diff is collapsed.
#ifndef Millepede_H
#define Millepede_H 1
#include "core/algorithm/Algorithm.h"
#include "objects/Track.h"
namespace corryvreckan {
/** @ingroup Algorithms
*
* Implementation of the Millepede algorithm.
*
* @author Christoph Hombach
* @date 2012-06-19
*/
class Millepede : public Algorithm {
public:
/// Constructor
Millepede(Configuration config, std::vector<Detector*> detectors);
/// Destructor
virtual ~Millepede();
void initialise();
void finalise();
StatusCode run(Clipboard*);
virtual void updateGeometry();
private:
struct Equation {
double rmeas;
double weight;
std::vector<int> indG;
std::vector<double> derG;
std::vector<int> indL;
std::vector<double> derL;
std::vector<int> indNL;
std::vector<double> derNL;
std::vector<double> slopes;
};
struct Constraint {
/// Right-hand side (Lagrange multiplier)
double rhs;
/// Coefficients
std::vector<double> coefficients;
};
/// (Re-)initialise matrices and vectors.
bool reset(const unsigned int nPlanes, const double startfact);
/// Setup the constraint equations.
void setConstraints(const unsigned int nPlanes);
/// Define a single constraint equation
void addConstraint(const std::vector<double>& dercs, const double rhs);
/// Add the equations for one track and do the local fit.
bool putTrack(Track* track, const unsigned int nPlanes);
/// Store the parameters for one measurement.
void addEquation(std::vector<Equation>& equations,
const std::vector<double>& derlc,
const std::vector<double>& dergb,
const std::vector<double>& dernl,
const std::vector<int>& dernli,
const std::vector<double>& dernls,
const double rmeas,
const double sigma);
// Perform local parameters fit using the equations for one track.
bool fitTrack(const std::vector<Equation>& equations,
std::vector<double>& trackParams,
const bool singlefit,
const unsigned int iteration);
// Perform global parameters fit.
bool fitGlobal();
/// Print the results of the global parameters fit.
bool printResults();
/// Matrix inversion and solution for global fit.
int invertMatrix(std::vector<std::vector<double>>& v, std::vector<double>& b, const int n);
// Matrix inversion and solution for local fit.
int invertMatrixLocal(std::vector<std::vector<double>>& v, std::vector<double>& b, const int n);
/// Return the limit in chi2 / ndof for n sigmas.
double chi2Limit(const int n, const int nd) const;
/// Multiply matrix and vector
bool multiplyAX(const std::vector<std::vector<double>>& a,
const std::vector<double>& x,
std::vector<double>& y,
const unsigned int n,
const unsigned int m);
/// Multiply matrices
bool multiplyAVAt(const std::vector<std::vector<double>>& v,
const std::vector<std::vector<double>>& a,
std::vector<std::vector<double>>& w,
const unsigned int n,
const unsigned int m);
Tracks m_alignmenttracks;
int m_numberOfTracksForAlignment;
/// Number of global derivatives
unsigned int m_nagb;
/// Number of local derivatives
unsigned int m_nalc = 4;
/// Degrees of freedom
std::vector<bool> m_dofs;
/// Default degrees of freedom
std::vector<bool> m_dofsDefault;
/// Equations for each track
std::vector<std::vector<Equation>> m_equations;
/// Constraint equations
std::vector<Constraint> m_constraints;
/// Flag for each global parameter whether it is fixed or not.
std::vector<bool> m_fixed;
/// Sigmas for each global parameter.
std::vector<double> m_psigm;
std::vector<std::vector<double>> m_cgmat;
std::vector<std::vector<double>> m_corrm;
std::vector<std::vector<double>> m_clcmat;
std::vector<double> m_bgvec;
std::vector<double> m_corrv;
std::vector<double> m_diag;
/// Difference in misalignment parameters with respect to initial values.
std::vector<double> m_dparm;
/// Mapping of internal numbering to geometry service planes.
std::map<std::string, unsigned int> m_millePlanes;
bool m_debug = false;
/// Flag to switch on/off iterations in the global fit.
bool m_iterate = true;
/// Residual cut after the first iteration.
double m_rescut;
/// Residual cut in the first iteration.
double m_rescut_init;
/// Factor in chisquare / ndof cut.
double m_cfactr;
/// Reference value for chisquare / ndof cut factor.
double m_cfactref;
// Number of standard deviations for chisquare / ndof cut.
int m_nstdev;
/// Number of "full" iterations (with geometry updates).
unsigned int m_nIterations;
/// Sigmas for each degree of freedom
std::vector<double> m_sigmas;
/// Planes to be kept fixed
std::vector<unsigned int> m_fixedPlanes;
/// Flag to fix all degrees of freedom or only the translations.
bool m_fix_all;
};
} // namespace corryvreckan
#endif // Millepede_H
## Millepede
**Maintainer**: Simon Spannagel (<simon.spannagel@cern.ch>)
**Status**: Work in progress
#### Description
This implementation of the Millepede algorithm has been taken from the [Kepler framework](https://gitlab.cern.ch/lhcb/Kepler) used for test beam data analysis within the LHCb collaboration. It has been written by Christoph Hombach and has seen contributions from Tim Evans and Heinrich Schindler. This version is only slightly adapted to the Corryvreckan framework.
The Millepede algorthm allows a simultaneous fit of both the tracks and the alignment constants.
#### Parameters
* `number_of_tracks`: Number of tracks used in the alignment method chosen. Default value is `20000`.
* `iterations`: Number of times the chosen alignment method is to be iterated. Default value is `3`.
* `dofs`: Degrees of freedom to be aligned. This parameter should be given as vector of six boolean values for the parameters "Translation X", "Translation Y", "Translation Z", "Rotation X", "Rotation Y" and "Rotation Z". The default setting is an alignment of all parameters except for "Translation Z", i.e. `dofs = true, true, false, true, true, true`.
* `residual_cut`: Residual cut to reject a track as an outlier. Default value is `0.05mm`;
* `residual_cut_init`: Initial residual cut for outlier rejection in the first iteration. This value is applied for the first iteration and replaced by `residual_cut` thereafter. Default value is `0.6mm`.
* `nstddev`: Cut to reject track candidates based on their Chi2/ndof value. Default value is `0`, i.e. the feature is disabled.
* `sigmas`: Uncertainties for each of the alignment parameters. Defaults to `0.05, 0.05, 0.5, 0.005, 0.005, 0.005`.
#### Plots produced
No plots are produced.
#### Usage
```toml
[Millepede]
iterations = 10
dofs = true, true, false, true, true, true
```
Parameters to be used in multiple algorithms can also be defined globally at the top of the configuration file. This is highly encouraged for parameters such as `DUT` and `reference`.
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment