Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DKinFitter.h
Go to the documentation of this file.
1 #ifndef _DKinFitter_
2 #define _DKinFitter_
3 
4 #include <algorithm>
5 #include <utility>
6 #include <math.h>
7 #include <iostream>
8 #include <map>
9 #include <set>
10 #include <limits>
11 
12 #include "TVector3.h"
13 #include "TMatrixD.h"
14 #include "TMatrixDSym.h"
15 #include "TMatrixFSym.h"
16 #include "TLorentzVector.h"
17 #include "TMath.h"
18 #include "TDecompLU.h"
19 
20 #include "DKinFitParticle.h"
21 #include "DKinFitConstraint.h"
22 #include "DKinFitConstraint_Mass.h"
23 #include "DKinFitConstraint_P4.h"
26 
27 using namespace std;
28 
30 {
36 };
37 
38 class DKinFitUtils; //forward declaration
39 
40 class DKinFitter //purely virtual: cannot directly instantiate class, can only inherit from it
41 {
42  public:
43 
44  /****************************************************************** PRIMARY *****************************************************************/
45 
46  //CONSTRUCTOR
47  DKinFitter(DKinFitUtils* locKinFitUtils);
48 
49  //RESET
50  void Reset_NewEvent(void);
51  void Reset_NewFit(void);
52 
53  //ADD CONSTRAINTS
54  void Add_Constraint(const shared_ptr<DKinFitConstraint>& locKinFitConstraint){dKinFitConstraints.insert(locKinFitConstraint);}
55  void Add_Constraints(const set<shared_ptr<DKinFitConstraint>>& locKinFitConstraints);
56 
57  //FIT!
58  bool Fit_Reaction(void); //IF YOU OVERRIDE THIS METHOD IN THE DERIVED CLASS, MAKE SURE YOU CALL THIS BASE CLASS METHOD!!
59 
60  //RECYCLE MEMORY //only call if you are discarding the results from the previous fit
61  void Recycle_LastFitMemory(void); //e.g. fit failed, or used to get a vertex guess
62 
63  /**************************************************************** FIT CONTROL ***************************************************************/
64 
65  //GET CONTROL VARIABLES
66  int Get_DebugLevel(void) const{return dDebugLevel;}
67  unsigned int Get_MaxNumIterations(void) const{return dMaxNumIterations;}
68  double Get_ConvergenceChiSqDiff(void) const{return dConvergenceChiSqDiff;}
69  double Get_ConvergenceChiSqDiff_LastResort(void) const{return dConvergenceChiSqDiff_LastResort;}
70 
71  //SET CONTROL VARIABLES
72  void Set_DebugLevel(int locDebugLevel);
73  void Set_MaxNumIterations(unsigned int locMaxNumIterations){dMaxNumIterations = locMaxNumIterations;}
74  void Set_ConvergenceChiSqDiff(double locConvergenceChiSqDiff){dConvergenceChiSqDiff = locConvergenceChiSqDiff;}
75  void Set_ConvergenceChiSqDiff_LastResort(double locConvergenceChiSqDiff){dConvergenceChiSqDiff_LastResort = locConvergenceChiSqDiff;}
76 
77  /************************************************************** GET FIT RESULTS *************************************************************/
78 
79  //GET STATUS (Fit_Reaction returns true/false for success/fail, but this provides more details on the failures
80  DKinFitStatus Get_KinFitStatus(void) const{return dKinFitStatus;}
81 
82  //GET FIT INFORMATION
83  unsigned int Get_NumUnknowns(void) const{return dNumXi;}
84  unsigned int Get_NumMeasurables(void) const{return dNumEta;}
85  unsigned int Get_NumConstraintEquations(void) const{return dNumF;}
86 
87  //GET FIT QUALITY RESULTS
88  double Get_ChiSq(void) const{return dChiSq;}
89  double Get_ConfidenceLevel(void) const{return dConfidenceLevel;}
90  unsigned int Get_NDF(void) const{return dNDF;}
91  void Get_Pulls(map<shared_ptr<DKinFitParticle>, map<DKinFitPullType, double> >& locPulls) const{locPulls = dPulls;} //key is particle, 2nd key is param type
92 
93  //GET UNCERTAINTIES
94  const TMatrixDSym& Get_VEta(void) {return dVEta;}
95  const TMatrixDSym& Get_VXi(void) {return dVXi;}
96  const TMatrixDSym& Get_V(void) {return dV;}
97 
98  //GET OUTPUT PARTICLES & CONSTRAINTS
99  set<shared_ptr<DKinFitConstraint>> Get_KinFitConstraints(void) const{return dKinFitConstraints;}
100  set<shared_ptr<DKinFitParticle>> Get_KinFitParticles(void) const{return dKinFitParticles;}
101 
102  /************************************************************ END GET FIT RESULTS ***********************************************************/
103 
104  private:
105 
106  //PRIVATE DEFAULT CONSTRUCTOR
107  DKinFitter(void){}; //Cannot use default constructor. Must construct with DKinFitUtils as argument
108 
109  /************************************************************ UTILITY FUNCTIONS *************************************************************/
110 
111  template <typename DType> set<shared_ptr<DType>> Get_Constraints(void) const;
112  template <typename DType> set<shared_ptr<DType>> Get_Constraints(const set<shared_ptr<DKinFitConstraint>>& locConstraints) const;
113 
114  //When we ask this question, we want to know, is the particle's information used in a constraint of the given type
115  //For example, if a decaying particle is used for a vertex constraint, the particle's that DEFINE the decaying particle are also included
116  template <typename DType> set<shared_ptr<DType>> Get_Constraints(const shared_ptr<DKinFitParticle>& locKinFitParticle, bool locOnlyDirectFlag = false) const;
117  template <typename DType> bool Get_IsInConstraint(const shared_ptr<DKinFitParticle>& locKinFitParticle, bool locOnlyDirectFlag = false) const;
118  template <typename DType> bool Get_IsIndirectlyInConstraint(const shared_ptr<DKinFitParticle>& locKinFitParticle) const;
119 
120  bool Get_IsConstrainingVertex(const shared_ptr<DKinFitParticle>& locKinFitParticle) const;
121  bool Get_IsTimeConstrained(const shared_ptr<DKinFitParticle>& locKinFitParticle) const;
122 
123  /*********************************************************** FIT INITIALIZATION *************************************************************/
124 
125  void Prepare_ConstraintsAndParticles(void);
126  void Set_MatrixSizes(void);
127  void Resize_Matrices(void);
128  void Zero_Matrices(void);
129  void Fill_InputMatrices(void);
130 
131  /************************************************************ CALCULATE MATRICES ************************************************************/
132 
133  bool Iterate(void);
134 
135  bool Calc_dS(void);
136  bool Calc_dU(void);
137  void Calc_dVdEta(void);
138 
139  void Calc_dF(void);
140 
141  void Calc_dF_P4(int locFIndex, const DKinFitParticle* locKinFitParticle, double locStateSignMultiplier);
142  void Calc_dF_MassDerivs(size_t locFIndex, const DKinFitParticle* locKinFitParticle, TLorentzVector locXP4, double locStateSignMultiplier, bool locIsConstrainedParticle);
143 
144  void Calc_dF_Vertex(size_t locFIndex, const DKinFitParticle* locKinFitParticle, const DKinFitParticle* locKinFitParticle_DecayingSource, double locStateSignMultiplier);
145  void Calc_dF_Vertex_NotDecaying(size_t locFIndex, const DKinFitParticle* locKinFitParticle);
146  void Calc_dF_Vertex_Decaying_Accel(size_t locFIndex, const DKinFitParticle* locKinFitParticle, const DKinFitParticle* locKinFitParticle_DecayingSource, double locStateSignMultiplier);
147  void Calc_dF_Vertex_Decaying_NonAccel(size_t locFIndex, const DKinFitParticle* locKinFitParticle, const DKinFitParticle* locKinFitParticle_DecayingSource, double locStateSignMultiplier);
148  void Calc_Vertex_Params(const DKinFitParticle* locKinFitParticle, double& locJ, TVector3& locQ, TVector3& locM, TVector3& locD);
149  TVector3 Calc_VertexParams_P4DerivedAtCommonVertex(const DKinFitParticle* locKinFitParticle);
150 
151  /************************************************************* UPDATE & FINAL ***************************************************************/
152 
153  void Update_ParticleParams(void);
154  void Calc_Pulls(void);
155  void Set_FinalTrackInfo(void);
156  void Update_CovarianceMatrices(bool locDecayingParticlesOnlyFlag);
157 
158  /***************************************************** FIT CONTROL AND UTILITY VARIABLES ****************************************************/
159 
161 
163  unsigned int dDebugLevel;
164 
165  unsigned int dMaxNumIterations;
166 
168  double dConvergenceChiSqDiff_LastResort; //if max # iterations hit, use this for final check (sometimes chisq walks (very slightly) forever without any meaningful change in the variables)
169 
170  /******************************************************** CONSTRAINTS AND PARTICLES *********************************************************/
171 
172  set<shared_ptr<DKinFitConstraint>> dKinFitConstraints;
173  set<shared_ptr<DKinFitParticle>> dKinFitParticles;
174  map<shared_ptr<DKinFitParticle>, set<shared_ptr<DKinFitConstraint>>> dParticleConstraintMap; //these particles are either directly or indirectly in these constraints
175  map<shared_ptr<DKinFitParticle>, set<shared_ptr<DKinFitConstraint>>> dParticleConstraintMap_Direct; //these particles are directly in these constraints
176  //indirect: if only present to define the momentum of a decaying particle
177 
178  /*********************************************************** FIT MATRIX VARIABLES ***********************************************************/
179 
180  unsigned int dNumXi; //num unknowns
181  unsigned int dNumEta; //num measurables
182  unsigned int dNumF; //num constraint eqs
183 
184  TMatrixD dXi; //unmeasurable unknowns
185  TMatrixD dEta; //observables
186  TMatrixD dY; //first approximation of observables (initial measurements)
187  TMatrixDSym dVY; //convariance matrix of dY
188 
189  TMatrixDSym dS; // the covariance matrix of the uncertainties of whether the constraint equations are satisfied
190  TMatrixDSym dS_Inverse;
191  TMatrixDSym dU;
192  TMatrixDSym dU_Inverse;
193 
194  TMatrixD dF; //constraint equations with eta dependence
195  TMatrixD dEpsilon; //dY - dEta
196  TMatrixD dLambda; //lagrange multipliers of constraint equations
197  TMatrixD dLambda_T;
198 
199  TMatrixD dF_dEta; //partial derivative of constraint equations wrst the observables
200  TMatrixD dF_dEta_T;
201  TMatrixD dF_dXi; //partial derivative of constraint equations wrst the unmeasurable unknowns
202  TMatrixD dF_dXi_T;
203 
204  TMatrixDSym dVXi; //covariance matrix of dXi
205  TMatrixDSym dVEta; //covariance matrix of dEta
206  TMatrixDSym dV; //full covariance matrix: dVEta at top-left and dVXi at bottom-right (+ the eta, xi covariance)
207 
208  /*************************************************************** FIT RESULTS ****************************************************************/
209 
210  double dChiSq;
211  unsigned int dNDF;
213 
214  //Pulls: 2nd dimension can be E, x, y, z, t for neutral showers; px, py, pz, x, y, z, t for charged or neutral tracks; or a subset of these
215  map<shared_ptr<DKinFitParticle>, map<DKinFitPullType, double> > dPulls; //key is particle, 2nd key is param type
216 };
217 
218 inline void DKinFitter::Add_Constraints(const set<shared_ptr<DKinFitConstraint>>& locKinFitConstraints)
219 {
220  dKinFitConstraints.insert(locKinFitConstraints.begin(), locKinFitConstraints.end());
221 }
222 
223 template <typename DType> inline set<shared_ptr<DType>> DKinFitter::Get_Constraints(void) const
224 {
225  //Get all constraints of a given type
226  return Get_Constraints<DType>(dKinFitConstraints);
227 }
228 
229 template <typename DType> inline set<shared_ptr<DType>> DKinFitter::Get_Constraints(const set<shared_ptr<DKinFitConstraint>>& locConstraints) const
230 {
231  //Get all constraints of a given type
232  set<shared_ptr<DType>> locTypeConstraints;
233  auto locConstraintIterator = locConstraints.begin();
234  for(; locConstraintIterator != locConstraints.end(); ++locConstraintIterator)
235  {
236  auto locConstraint = std::dynamic_pointer_cast<DType>(*locConstraintIterator);
237  if(locConstraint != NULL)
238  locTypeConstraints.insert(locConstraint);
239  }
240  return locTypeConstraints;
241 }
242 
243 template <typename DType> inline bool DKinFitter::Get_IsInConstraint(const shared_ptr<DKinFitParticle>& locKinFitParticle, bool locOnlyDirectFlag) const
244 {
245  //Return whether a particle is in a constraint of the given type
246  return !Get_Constraints<DType>(locKinFitParticle, locOnlyDirectFlag).empty();
247 }
248 
249 template <typename DType> inline bool DKinFitter::Get_IsIndirectlyInConstraint(const shared_ptr<DKinFitParticle>& locKinFitParticle) const
250 {
251  //Return whether a particle is indirectly in a constraint of the given type //in it, but not directly
252  return (Get_IsInConstraint<DType>(locKinFitParticle, false) && !Get_IsInConstraint<DType>(locKinFitParticle, true));
253 }
254 
255 template <typename DType> inline set<shared_ptr<DType>> DKinFitter::Get_Constraints(const shared_ptr<DKinFitParticle>& locKinFitParticle, bool locOnlyDirectFlag) const
256 {
257  //Get all constraints of a given type that a given particle is in: either directly or INDIRECTLY
258  //If the particle is used to define the p3 of a decaying particle that is used in a given constraint, it is included
259  auto& locConstraintMap = locOnlyDirectFlag ? dParticleConstraintMap_Direct : dParticleConstraintMap;
260  auto locMapIterator = locConstraintMap.find(locKinFitParticle);
261  if(locMapIterator == locConstraintMap.end())
262  return set<shared_ptr<DType>>();
263 
264  auto& locParticleConstraints = locMapIterator->second;
265  return Get_Constraints<DType>(locParticleConstraints);
266 }
267 
268 #endif // _DKinFitter_
TMatrixDSym dS_Inverse
Definition: DKinFitter.h:190
TMatrixD dXi
Definition: DKinFitter.h:184
TMatrixD dEta
Definition: DKinFitter.h:185
DKinFitStatus
Definition: DKinFitter.h:29
TMatrixD dF
Definition: DKinFitter.h:194
map< shared_ptr< DKinFitParticle >, set< shared_ptr< DKinFitConstraint > > > dParticleConstraintMap
Definition: DKinFitter.h:174
unsigned int dMaxNumIterations
Definition: DKinFitter.h:165
double dConvergenceChiSqDiff_LastResort
Definition: DKinFitter.h:168
double dConvergenceChiSqDiff
Definition: DKinFitter.h:167
unsigned int dDebugLevel
Definition: DKinFitter.h:163
double dChiSq
Definition: DKinFitter.h:210
TMatrixDSym dVXi
Definition: DKinFitter.h:204
unsigned int dNumF
Definition: DKinFitter.h:182
unsigned int dNumXi
Definition: DKinFitter.h:180
TMatrixDSym dS
Definition: DKinFitter.h:189
TMatrixD dLambda_T
Definition: DKinFitter.h:197
void Add_Constraint(const shared_ptr< DKinFitConstraint > &locKinFitConstraint)
Definition: DKinFitter.h:54
void Set_ConvergenceChiSqDiff_LastResort(double locConvergenceChiSqDiff)
Definition: DKinFitter.h:75
const TMatrixDSym & Get_V(void)
Definition: DKinFitter.h:96
unsigned int dNDF
Definition: DKinFitter.h:211
TMatrixDSym dU_Inverse
Definition: DKinFitter.h:192
bool Get_IsIndirectlyInConstraint(const shared_ptr< DKinFitParticle > &locKinFitParticle) const
Definition: DKinFitter.h:249
TMatrixD dY
Definition: DKinFitter.h:186
TMatrixD dLambda
Definition: DKinFitter.h:196
unsigned int Get_NumMeasurables(void) const
Definition: DKinFitter.h:84
TMatrixDSym dVEta
Definition: DKinFitter.h:205
DKinFitter(void)
Definition: DKinFitter.h:107
double dConfidenceLevel
Definition: DKinFitter.h:212
unsigned int Get_NumConstraintEquations(void) const
Definition: DKinFitter.h:85
bool Get_IsInConstraint(const shared_ptr< DKinFitParticle > &locKinFitParticle, bool locOnlyDirectFlag=false) const
Definition: DKinFitter.h:243
TMatrixD dF_dEta
Definition: DKinFitter.h:199
map< shared_ptr< DKinFitParticle >, map< DKinFitPullType, double > > dPulls
Definition: DKinFitter.h:215
void Set_MaxNumIterations(unsigned int locMaxNumIterations)
Definition: DKinFitter.h:73
const TMatrixDSym & Get_VEta(void)
Definition: DKinFitter.h:94
set< shared_ptr< DKinFitConstraint > > Get_KinFitConstraints(void) const
Definition: DKinFitter.h:99
unsigned int Get_NDF(void) const
Definition: DKinFitter.h:90
TMatrixDSym dVY
Definition: DKinFitter.h:187
double Get_ChiSq(void) const
Definition: DKinFitter.h:88
DKinFitUtils * dKinFitUtils
Definition: DKinFitter.h:160
TMatrixD dEpsilon
Definition: DKinFitter.h:195
TMatrixD dF_dXi
Definition: DKinFitter.h:201
void Add_Constraints(const set< shared_ptr< DKinFitConstraint >> &locKinFitConstraints)
Definition: DKinFitter.h:218
void Get_Pulls(map< shared_ptr< DKinFitParticle >, map< DKinFitPullType, double > > &locPulls) const
Definition: DKinFitter.h:91
TMatrixDSym dU
Definition: DKinFitter.h:191
const TMatrixDSym & Get_VXi(void)
Definition: DKinFitter.h:95
double Get_ConvergenceChiSqDiff_LastResort(void) const
Definition: DKinFitter.h:69
set< shared_ptr< DKinFitParticle > > Get_KinFitParticles(void) const
Definition: DKinFitter.h:100
double Get_ConfidenceLevel(void) const
Definition: DKinFitter.h:89
DKinFitStatus Get_KinFitStatus(void) const
Definition: DKinFitter.h:80
int Get_DebugLevel(void) const
Definition: DKinFitter.h:66
TMatrixD dF_dXi_T
Definition: DKinFitter.h:202
void Set_ConvergenceChiSqDiff(double locConvergenceChiSqDiff)
Definition: DKinFitter.h:74
set< shared_ptr< DKinFitConstraint > > dKinFitConstraints
Definition: DKinFitter.h:172
unsigned int Get_MaxNumIterations(void) const
Definition: DKinFitter.h:67
unsigned int Get_NumUnknowns(void) const
Definition: DKinFitter.h:83
set< shared_ptr< DKinFitParticle > > dKinFitParticles
Definition: DKinFitter.h:173
DKinFitStatus dKinFitStatus
Definition: DKinFitter.h:162
unsigned int dNumEta
Definition: DKinFitter.h:181
TMatrixD dF_dEta_T
Definition: DKinFitter.h:200
set< shared_ptr< DType > > Get_Constraints(void) const
Definition: DKinFitter.h:223
map< shared_ptr< DKinFitParticle >, set< shared_ptr< DKinFitConstraint > > > dParticleConstraintMap_Direct
Definition: DKinFitter.h:175
double Get_ConvergenceChiSqDiff(void) const
Definition: DKinFitter.h:68
TMatrixDSym dV
Definition: DKinFitter.h:206