GblTrajectory.h 5.68 KB
Newer Older
1
2
3
4
5
6
7
/*
 * GblTrajectory.h
 *
 *  Created on: Aug 18, 2011
 *      Author: kleinwrt
 */

8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
/** \file
 *  GblTrajectory definition.
 *
 *  \author Claus Kleinwort, DESY, 2011 (Claus.Kleinwort@desy.de)
 *
 *  \copyright
 *  Copyright (c) 2011 - 2016 Deutsches Elektronen-Synchroton,
 *  Member of the Helmholtz Association, (DESY), HAMBURG, GERMANY \n\n
 *  This library is free software; you can redistribute it and/or modify
 *  it under the terms of the GNU Library General Public License as
 *  published by the Free Software Foundation; either version 2 of the
 *  License, or (at your option) any later version. \n\n
 *  This library 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 Library General Public License for more details. \n\n
 *  You should have received a copy of the GNU Library General Public
 *  License along with this program (see the file COPYING.LIB for more
 *  details); if not, write to the Free Software Foundation, Inc.,
 *  675 Mass Ave, Cambridge, MA 02139, USA.
 */

30
31
32
33
34
#ifndef GBLTRAJECTORY_H_
#define GBLTRAJECTORY_H_

#include "GblPoint.h"
#include "GblData.h"
35
#include "GblPoint.h"
36
37
38
39
#include "BorderedBandMatrix.h"
#include "MilleBinary.h"
#include "TMatrixDSymEigen.h"

40
//! Namespace for the general broken lines package
Claus Kleinwort's avatar
Claus Kleinwort committed
41
namespace gbl {
42

43
44
45
46
47
48
49
/// GBL trajectory.
/**
 * List of GblPoints ordered by arc length.
 * Can be fitted and optionally written to MP-II binary file.
 */
class GblTrajectory {
public:
50
	GblTrajectory(const std::vector<GblPoint> &aPointList, bool flagCurv = true,
Claus Kleinwort's avatar
Claus Kleinwort committed
51
			bool flagU1dir = true, bool flagU2dir = true);
52
	GblTrajectory(const std::vector<GblPoint> &aPointList, unsigned int aLabel,
53
54
			const TMatrixDSym &aSeed, bool flagCurv = true, bool flagU1dir =
					true, bool flagU2dir = true);
55
56
57
58
59
60
	GblTrajectory(
			const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointaAndTransList);
	GblTrajectory(
			const std::vector<std::pair<std::vector<GblPoint>, TMatrixD> > &aPointaAndTransList,
			const TMatrixD &extDerivatives, const TVectorD &extMeasurements,
			const TVectorD &extPrecisions);
61
	virtual ~GblTrajectory();
62
	bool isValid() const;
63
	unsigned int getNumPoints() const;
64
	unsigned int getResults(int aSignedLabel, TVectorD &localPar,
65
			TMatrixDSym &localCov) const;
66
67
68
69
70
71
	unsigned int getMeasResults(unsigned int aLabel, unsigned int &numRes,
			TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors,
			TVectorD &aDownWeights);
	unsigned int getScatResults(unsigned int aLabel, unsigned int &numRes,
			TVectorD &aResiduals, TVectorD &aMeasErrors, TVectorD &aResErrors,
			TVectorD &aDownWeights);
Claus Kleinwort's avatar
Claus Kleinwort committed
72
73
	void getLabels(std::vector<unsigned int> &aLabelList);
	void getLabels(std::vector<std::vector< unsigned int> > &aLabelList);
74
	unsigned int fit(double &Chi2, int &Ndf, double &lostWeight,
75
76
			std::string optionList = "");
	void milleOut(MilleBinary &aMille);
77
78
79
	void printTrajectory(unsigned int level = 0);
	void printPoints(unsigned int level = 0);
	void printData();
80
81

private:
82
83
84
	unsigned int numAllPoints; ///< Number of all points on trajectory
	std::vector<unsigned int> numPoints; ///< Number of points on (sub)trajectory
	unsigned int numTrajectories; ///< Number of trajectories (in composed trajectory)
85
	unsigned int numOffsets; ///< Number of (points with) offsets on trajectory
86
87
	unsigned int numInnerTrans; ///< Number of inner transformations to external parameters
	unsigned int numCurvature; ///< Number of curvature parameters (0 or 1) or external parameters
88
89
	unsigned int numParameters; ///< Number of fit parameters
	unsigned int numLocals; ///< Total number of (additional) local parameters
90
	unsigned int numMeasurements; ///< Total number of measurements
91
	unsigned int externalPoint; ///< Label of external point (or 0)
92
	bool constructOK; ///< Trajectory has been successfully constructed (ready for fit/output)
93
	bool fitOK; ///< Trajectory has been successfully fitted (results are valid)
94
	std::vector<unsigned int> theDimension; ///< List of active dimensions (0=u1, 1=u2) in fit
95
	std::vector<std::vector<GblPoint> > thePoints; ///< (list of) List of points on trajectory
96
	std::vector<GblData> theData; ///< List of data blocks
97
98
	std::vector<unsigned int> measDataIndex; ///< mapping points to data blocks from measurements
	std::vector<unsigned int> scatDataIndex; ///< mapping points to data blocks from scatterers
99
	TMatrixDSym externalSeed; ///< Precision (inverse covariance matrix) of external seed
100
101
102
103
104
	std::vector<TMatrixD> innerTransformations; ///< Transformations at innermost points of
	// composed trajectory (from common external parameters)
	TMatrixD externalDerivatives; // Derivatives for external measurements of composed trajectory
	TVectorD externalMeasurements; // Residuals for external measurements of composed trajectory
	TVectorD externalPrecisions; // Precisions for external measurements of composed trajectory
105
	VVector theVector; ///< Vector of linear equation system
106
107
108
	BorderedBandMatrix theMatrix; ///< (Bordered band) matrix of linear equation system

	std::pair<std::vector<unsigned int>, TMatrixD> getJacobian(
109
			int aSignedLabel) const;
110
	void getFitToLocalJacobian(std::vector<unsigned int> &anIndex,
111
			SMatrix55 &aJacobian, const GblPoint &aPoint, unsigned int measDim,
112
			unsigned int nJacobian = 1) const;
113
	void getFitToKinkJacobian(std::vector<unsigned int> &anIndex,
114
			SMatrix27 &aJacobian, const GblPoint &aPoint) const;
115
	void construct();
116
117
118
	void defineOffsets();
	void calcJacobians();
	void prepare();
119
	void buildLinearEquationSystem();
120
121
	void predict();
	double downWeight(unsigned int aMethod);
122
123
	void getResAndErr(unsigned int aData, double &aResidual,
			double &aMeadsError, double &aResError, double &aDownWeight);
124
};
125
}
126
#endif /* GBLTRAJECTORY_H_ */