CMS 3D CMS Logo

List of all members | Public Types | Public Member Functions | Protected Member Functions | Protected Attributes | Private Member Functions | Private Attributes
tauImpactParameter::LagrangeMultipliersFitter Class Referenceabstract

#include <LagrangeMultipliersFitter.h>

Public Types

enum  ConvergeProc { ConstraintMin = 0, Chi2Min, Chi2AndConstaintMin }
 
enum  Parameters {
  par_vx = 0, par_vy, par_vz, par_px,
  par_py, par_pz, par_m, npardim
}
 
enum  Position { pos_x = 0, pos_y, pos_z, nposdim }
 

Public Member Functions

virtual double chiSquare ()
 
virtual double cSum ()
 
virtual bool fit ()
 
virtual LorentzVectorParticle getMother ()=0
 
virtual std::vector< LorentzVectorParticlegetRefitDaughters ()=0
 
virtual bool isConfigured ()
 
virtual bool isConverged ()
 
 LagrangeMultipliersFitter ()
 
virtual double nConstraints ()=0
 
virtual int nDaughters ()=0
 
virtual double ndf ()=0
 
virtual double nIter ()
 
virtual void setMaxDelta (double MaxDelta)
 
virtual void setNIterMax (int Nitermax)
 
virtual void setWeight (double weight)
 
virtual ~LagrangeMultipliersFitter ()
 

Protected Member Functions

virtual TVectorD value (const TVectorD &v)=0
 

Protected Attributes

TMatrixTSym< double > cov_
 
TMatrixTSym< double > cov_0_
 
bool isConfigured_
 
bool isFit_
 
TVectorD par_
 
TVectorD par_0_
 

Private Member Functions

bool applyLagrangianConstraints ()
 
double chiSquare (const TVectorT< double > &delta_alpha, const TVectorT< double > &lambda, const TMatrixT< double > &D, const TVectorT< double > &d)
 
double chiSquareUsingInitalPoint (const TVectorT< double > &alpha, const TVectorT< double > &lambda)
 
TMatrixT< double > computeVariance ()
 
double constraintDelta (const TVectorT< double > &par)
 
TMatrixT< double > derivative ()
 

Private Attributes

double chi2_
 
double chi2prev_
 
TMatrixT< double > D_
 
double delta_
 
double epsilon_
 
double maxDelta_
 
double niter_
 
double nitermax_
 
double ScaleFactor_
 
TMatrixTSym< double > V_alpha0_inv_
 
TMatrixT< double > V_corr_prev_
 
TMatrixTSym< double > V_D_
 
double weight_
 

Detailed Description

Definition at line 18 of file LagrangeMultipliersFitter.h.

Member Enumeration Documentation

◆ ConvergeProc

◆ Parameters

◆ Position

Constructor & Destructor Documentation

◆ LagrangeMultipliersFitter()

LagrangeMultipliersFitter::LagrangeMultipliersFitter ( )

◆ ~LagrangeMultipliersFitter()

virtual tauImpactParameter::LagrangeMultipliersFitter::~LagrangeMultipliersFitter ( )
inlinevirtual

Definition at line 25 of file LagrangeMultipliersFitter.h.

25 {};

Member Function Documentation

◆ applyLagrangianConstraints()

bool LagrangeMultipliersFitter::applyLagrangianConstraints ( )
private

Definition at line 50 of file LagrangeMultipliersFitter.cc.

References alpha, correctionTermsCaloMet_cff::C, chi2_, Chi2AndConstaintMin, chi2prev_, chiSquareUsingInitalPoint(), constraintDelta(), ConstraintMin, cov_0_, ztail::d, D_, delta_, derivative(), change_name::diff, GeomDetEnumerators::DT, epsilon_, createfilelist::int, MainPageGenerator::l, maxDelta_, nConstraints(), par_, par_0_, alignCSCRings::s, ScaleFactor_, V_D_, and value().

Referenced by fit().

50  {
51  if (V_D_.GetNrows() != nConstraints())
52  V_D_.ResizeTo(nConstraints(), nConstraints());
53  if (D_.GetNrows() != nConstraints() || D_.GetNcols() != par_.GetNrows())
54  D_.ResizeTo(nConstraints(), par_.GetNrows());
55 
56  // Setup intial values
57  TVectorT<double> alpha_A = par_;
58  TVectorT<double> alpha_0 = par_0_;
59  TVectorT<double> delta_alpha_A = alpha_A - alpha_0;
60  D_ = derivative();
61  TVectorT<double> d = value(par_);
62  TVectorT<double> C = D_ * delta_alpha_A - d;
63  TMatrixTSym<double> V_alpha0 = cov_0_;
64  TMatrixTSym<double> V_D_inv = V_alpha0;
65  V_D_inv.Similarity(D_);
66  double det = V_D_inv.Determinant();
67  TDecompBK Inverter(V_D_inv);
68  if (fabs(det) > 1e40) {
69  edm::LogWarning("LagrangeMultipliersFitter::Fit")
70  << "Fit failed: unable to invert SYM gain matrix LARGE Determinant" << det << " \n"
71  << std::endl;
72  return false;
73  }
74  if (!Inverter.Decompose()) {
75  edm::LogWarning("LagrangeMultipliersFitter::Fit") << "Fit failed: unable to invert SYM gain matrix " << det << " \n"
76  << std::endl;
77  return false;
78  }
79  V_D_ = Inverter.Invert();
80 
81  // solve equations
82  TVectorT<double> lambda = -1.0 * V_D_ * C;
83  TMatrixT<double> DT = D_;
84  DT.T();
85  TVectorT<double> alpha = alpha_0 - V_alpha0 * DT * lambda;
86 
87  // do while loop to see if the convergance criteria are satisfied
88  double s(1), stepscale(0.01);
89  chi2prev_ = chi2_;
90  double currentchi2(chiSquareUsingInitalPoint(alpha_A, lambda)), currentdelta(constraintDelta(par_));
91  TVectorT<double> alpha_s = alpha;
92  // convergence in 2 step procedure to minimize chi2 within MaxDelta_ of the constriants
93  // 1) Get within 5x MaxDelta_
94  // 2) converge based on improving chi2 and constrianed delta
95  unsigned int Proc = ConstraintMin;
96  if (constraintDelta(par_) < 5 * maxDelta_)
97  Proc = Chi2AndConstaintMin;
98  int NIter = (int)(1.0 / stepscale);
99  for (int iter = 0; iter < NIter; iter++) {
100  // compute safty cutoff for numberical constraint
101  double diff = 0;
102  for (int l = 0; l < alpha_s.GetNrows(); l++) {
103  if (diff < alpha_s(l) - alpha_A(l))
104  diff = alpha_s(l) - alpha_A(l);
105  }
106  double delta_alpha_s = constraintDelta(alpha_s);
107  if (Proc == ConstraintMin) {
108  if (delta_alpha_s < currentdelta || iter == NIter || diff < 100 * epsilon_) {
109  currentchi2 = chiSquareUsingInitalPoint(alpha_s, lambda);
110  currentdelta = delta_alpha_s;
111  ScaleFactor_ = s;
112  break;
113  }
114  } else if (Proc == Chi2AndConstaintMin) {
115  double chi2_s = chiSquareUsingInitalPoint(alpha_s, lambda);
116  if ((delta_alpha_s < currentdelta /*+maxDelta_*/ && chi2_s < currentchi2) || iter == NIter ||
117  diff < 100 * epsilon_) {
118  currentchi2 = chi2_s;
119  currentdelta = delta_alpha_s;
120  ScaleFactor_ = s;
121  break;
122  }
123  }
124  s -= stepscale;
125  alpha_s = alpha_A + s * (alpha - alpha_A);
126  }
127  // set chi2
128  chi2_ = currentchi2;
129  //set delta
130  delta_ = currentdelta;
131  par_ = alpha_s;
132  return true;
133 }
float alpha
Definition: AMPTWrapper.h:105
virtual TVectorD value(const TVectorD &v)=0
d
Definition: ztail.py:151
double chiSquareUsingInitalPoint(const TVectorT< double > &alpha, const TVectorT< double > &lambda)
double constraintDelta(const TVectorT< double > &par)
Log< level::Warning, false > LogWarning

◆ chiSquare() [1/2]

virtual double tauImpactParameter::LagrangeMultipliersFitter::chiSquare ( )
inlinevirtual

Definition at line 34 of file LagrangeMultipliersFitter.h.

References chi2_.

◆ chiSquare() [2/2]

double LagrangeMultipliersFitter::chiSquare ( const TVectorT< double > &  delta_alpha,
const TVectorT< double > &  lambda,
const TMatrixT< double > &  D,
const TVectorT< double > &  d 
)
private

Definition at line 162 of file LagrangeMultipliersFitter.cc.

References ztail::d.

165  {
166  double c2 = lambda * (D * delta_alpha + d);
167  return c2;
168 }
d
Definition: ztail.py:151
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:141

◆ chiSquareUsingInitalPoint()

double LagrangeMultipliersFitter::chiSquareUsingInitalPoint ( const TVectorT< double > &  alpha,
const TVectorT< double > &  lambda 
)
private

Definition at line 170 of file LagrangeMultipliersFitter.cc.

References alpha, cov_0_, mps_fire::i, dqmiolumiharvest::j, par_, par_0_, V_alpha0_inv_, and value().

Referenced by applyLagrangianConstraints().

171  {
172  if (cov_0_.GetNrows() != V_alpha0_inv_.GetNrows()) {
173  TMatrixTSym<double> V_alpha0 = cov_0_;
174  V_alpha0_inv_.ResizeTo(cov_0_.GetNrows(), cov_0_.GetNrows());
175  TDecompBK Inverter(V_alpha0);
176  if (!Inverter.Decompose()) { // handle rare case where inversion is not possible (ie assume diagonal)
177  edm::LogWarning("LagrangeMultipliersFitter::chiSquareUsingInitalPoint")
178  << "Error non-invertable Matrix... Calculating under assumption that correlations can be neglected!!!"
179  << std::endl;
180  for (int j = 0; j < par_.GetNrows(); j++) {
181  for (int i = 0; i < par_.GetNrows(); i++) {
182  if (i == j)
183  V_alpha0_inv_(i, j) = 1.0 / V_alpha0(i, j);
184  else
185  V_alpha0_inv_(i, j) = 0.0;
186  }
187  }
188  } else {
189  V_alpha0_inv_ = Inverter.Invert();
190  }
191  }
192 
193  TVectorT<double> alpha_0 = par_0_;
194  TVectorT<double> dalpha = alpha - alpha_0;
195  double c2_var = dalpha * (V_alpha0_inv_ * dalpha);
196  const TVectorT<double>& alpha_v = alpha;
197  double c2_constraints = lambda * value(alpha_v);
198  double c2 = c2_var + c2_constraints;
199  return c2;
200 }
float alpha
Definition: AMPTWrapper.h:105
virtual TVectorD value(const TVectorD &v)=0
Log< level::Warning, false > LogWarning

◆ computeVariance()

TMatrixT< double > LagrangeMultipliersFitter::computeVariance ( )
private

Definition at line 211 of file LagrangeMultipliersFitter.cc.

References cov_, cov_0_, D_, mps_fire::i, dqmiolumiharvest::j, ScaleFactor_, V_corr_prev_, and V_D_.

Referenced by fit().

211  {
212  TMatrixTSym<double> V_alpha0 = cov_0_;
213  TMatrixTSym<double> DTV_DD = V_D_.SimilarityT(D_);
214  TMatrixT<double> DTV_DDV = DTV_DD * V_alpha0;
215  TMatrixT<double> VDTV_DDV = V_alpha0 * DTV_DDV;
216  TMatrixT<double> CovCor = VDTV_DDV;
217  //CovCor*=ScaleFactor_;
218  if (V_corr_prev_.GetNrows() != V_alpha0.GetNrows()) {
219  V_corr_prev_.ResizeTo(V_alpha0.GetNrows(), V_alpha0.GetNrows());
220  V_corr_prev_ = CovCor;
221  } else {
222  V_corr_prev_ *= (1 - ScaleFactor_);
223  CovCor += V_corr_prev_;
224  V_corr_prev_ = CovCor;
225  }
226 
227  TMatrixT<double> V_alpha = V_alpha0 - CovCor;
228  for (int i = 0; i < cov_.GetNrows(); i++) {
229  for (int j = 0; j <= i; j++) {
230  cov_(i, j) = V_alpha(i, j);
231  }
232  }
233  return cov_;
234 }

◆ constraintDelta()

double LagrangeMultipliersFitter::constraintDelta ( const TVectorT< double > &  par)
private

Definition at line 202 of file LagrangeMultipliersFitter.cc.

References mps_fire::i, and value().

Referenced by applyLagrangianConstraints().

202  {
203  TVectorD d_par = value(par);
204  double delta_d(0);
205  for (int i = 0; i < d_par.GetNrows(); i++) {
206  delta_d += fabs(d_par(i));
207  }
208  return delta_d;
209 }
virtual TVectorD value(const TVectorD &v)=0

◆ cSum()

virtual double tauImpactParameter::LagrangeMultipliersFitter::cSum ( )
inlinevirtual

Definition at line 35 of file LagrangeMultipliersFitter.h.

References delta_.

◆ derivative()

TMatrixD LagrangeMultipliersFitter::derivative ( )
private

Definition at line 135 of file LagrangeMultipliersFitter.cc.

References epsilon_, mps_fire::i, dqmiolumiharvest::j, nConstraints(), par_, and value().

Referenced by applyLagrangianConstraints().

135  { // alway evaluated at current par
136  TMatrixD Derivatives(nConstraints(), par_.GetNrows());
137  TVectorD par_plus(par_.GetNrows());
138  TVectorD value_par(nConstraints());
139  TVectorD value_par_plus(nConstraints());
140  for (int j = 0; j < par_.GetNrows(); j++) {
141  for (int i = 0; i < par_.GetNrows(); i++) {
142  par_plus(i) = par_(i);
143  if (i == j)
144  par_plus(i) = par_(i) + epsilon_;
145  }
146  value_par = value(par_);
147  value_par_plus = value(par_plus);
148  for (int i = 0; i < nConstraints(); i++) {
149  Derivatives(i, j) = (value_par_plus(i) - value_par(i)) / epsilon_;
150  }
151  }
152  return Derivatives;
153 }
virtual TVectorD value(const TVectorD &v)=0
AlgebraicMatrix Derivatives
Definition: Definitions.h:35

◆ fit()

bool LagrangeMultipliersFitter::fit ( void  )
virtual

Definition at line 24 of file LagrangeMultipliersFitter.cc.

References applyLagrangianConstraints(), computeVariance(), cov_, cov_0_, delta_, isConfigured_, isConverged(), isFit_, maxDelta_, niter_, nitermax_, par_0_, and TriggerAnalyzer::passed.

Referenced by trackingPlots.Iteration::modules().

24  {
25  if (cov_.GetNrows() != par_0_.GetNrows()) {
26  // set cov to cov_0 until value is computed
27  cov_.ResizeTo(par_0_.GetNrows(), par_0_.GetNrows());
28  cov_ = cov_0_;
29  }
30  if (!isConfigured_)
31  return false;
32  if (isFit_)
33  return isConverged();
34  isFit_ = true;
35  niter_ = 0;
36  for (niter_ = 0; niter_ <= nitermax_; niter_++) {
38  if (!passed || (niter_ == nitermax_ && delta_ >= 4.0 * maxDelta_)) {
39  edm::LogWarning("LagrangeMultipliersFitter::Fit")
40  << "Reached Maximum number of iterations..." << niter_ << std::endl;
41  return false;
42  }
43  if (isConverged())
44  break;
45  }
47  return true;
48 }
Log< level::Warning, false > LogWarning

◆ getMother()

virtual LorentzVectorParticle tauImpactParameter::LagrangeMultipliersFitter::getMother ( )
pure virtual

◆ getRefitDaughters()

virtual std::vector<LorentzVectorParticle> tauImpactParameter::LagrangeMultipliersFitter::getRefitDaughters ( )
pure virtual

◆ isConfigured()

virtual bool tauImpactParameter::LagrangeMultipliersFitter::isConfigured ( )
inlinevirtual

◆ isConverged()

bool LagrangeMultipliersFitter::isConverged ( )
virtual

Definition at line 155 of file LagrangeMultipliersFitter.cc.

References delta_, and maxDelta_.

Referenced by fit().

155  {
156  if (delta_ < maxDelta_) {
157  return true;
158  }
159  return false;
160 }

◆ nConstraints()

virtual double tauImpactParameter::LagrangeMultipliersFitter::nConstraints ( )
pure virtual

◆ nDaughters()

virtual int tauImpactParameter::LagrangeMultipliersFitter::nDaughters ( )
pure virtual

◆ ndf()

virtual double tauImpactParameter::LagrangeMultipliersFitter::ndf ( )
pure virtual

◆ nIter()

virtual double tauImpactParameter::LagrangeMultipliersFitter::nIter ( )
inlinevirtual

Definition at line 36 of file LagrangeMultipliersFitter.h.

References niter_.

◆ setMaxDelta()

virtual void tauImpactParameter::LagrangeMultipliersFitter::setMaxDelta ( double  MaxDelta)
inlinevirtual

Definition at line 28 of file LagrangeMultipliersFitter.h.

References maxDelta_.

◆ setNIterMax()

virtual void tauImpactParameter::LagrangeMultipliersFitter::setNIterMax ( int  Nitermax)
inlinevirtual

Definition at line 29 of file LagrangeMultipliersFitter.h.

References nitermax_.

◆ setWeight()

virtual void tauImpactParameter::LagrangeMultipliersFitter::setWeight ( double  weight)
inlinevirtual

◆ value()

virtual TVectorD tauImpactParameter::LagrangeMultipliersFitter::value ( const TVectorD &  v)
protectedpure virtual

Member Data Documentation

◆ chi2_

double tauImpactParameter::LagrangeMultipliersFitter::chi2_
private

Definition at line 69 of file LagrangeMultipliersFitter.h.

Referenced by applyLagrangianConstraints(), and chiSquare().

◆ chi2prev_

double tauImpactParameter::LagrangeMultipliersFitter::chi2prev_
private

Definition at line 69 of file LagrangeMultipliersFitter.h.

Referenced by applyLagrangianConstraints().

◆ cov_

TMatrixTSym<double> tauImpactParameter::LagrangeMultipliersFitter::cov_
protected

Definition at line 50 of file LagrangeMultipliersFitter.h.

Referenced by computeVariance(), and fit().

◆ cov_0_

TMatrixTSym<double> tauImpactParameter::LagrangeMultipliersFitter::cov_0_
protected

◆ D_

TMatrixT<double> tauImpactParameter::LagrangeMultipliersFitter::D_
private

Definition at line 73 of file LagrangeMultipliersFitter.h.

Referenced by applyLagrangianConstraints(), and computeVariance().

◆ delta_

double tauImpactParameter::LagrangeMultipliersFitter::delta_
private

Definition at line 69 of file LagrangeMultipliersFitter.h.

Referenced by applyLagrangianConstraints(), cSum(), fit(), and isConverged().

◆ epsilon_

double tauImpactParameter::LagrangeMultipliersFitter::epsilon_
private

Definition at line 66 of file LagrangeMultipliersFitter.h.

Referenced by applyLagrangianConstraints(), and derivative().

◆ isConfigured_

bool tauImpactParameter::LagrangeMultipliersFitter::isConfigured_
protected

Definition at line 51 of file LagrangeMultipliersFitter.h.

Referenced by fit(), and isConfigured().

◆ isFit_

bool tauImpactParameter::LagrangeMultipliersFitter::isFit_
protected

Definition at line 52 of file LagrangeMultipliersFitter.h.

Referenced by fit().

◆ maxDelta_

double tauImpactParameter::LagrangeMultipliersFitter::maxDelta_
private

◆ niter_

double tauImpactParameter::LagrangeMultipliersFitter::niter_
private

Definition at line 69 of file LagrangeMultipliersFitter.h.

Referenced by fit(), and nIter().

◆ nitermax_

double tauImpactParameter::LagrangeMultipliersFitter::nitermax_
private

Definition at line 66 of file LagrangeMultipliersFitter.h.

Referenced by fit(), and setNIterMax().

◆ par_

TVectorD tauImpactParameter::LagrangeMultipliersFitter::par_
protected

◆ par_0_

TVectorD tauImpactParameter::LagrangeMultipliersFitter::par_0_
protected

◆ ScaleFactor_

double tauImpactParameter::LagrangeMultipliersFitter::ScaleFactor_
private

Definition at line 75 of file LagrangeMultipliersFitter.h.

Referenced by applyLagrangianConstraints(), and computeVariance().

◆ V_alpha0_inv_

TMatrixTSym<double> tauImpactParameter::LagrangeMultipliersFitter::V_alpha0_inv_
private

Definition at line 72 of file LagrangeMultipliersFitter.h.

Referenced by chiSquareUsingInitalPoint().

◆ V_corr_prev_

TMatrixT<double> tauImpactParameter::LagrangeMultipliersFitter::V_corr_prev_
private

Definition at line 76 of file LagrangeMultipliersFitter.h.

Referenced by computeVariance().

◆ V_D_

TMatrixTSym<double> tauImpactParameter::LagrangeMultipliersFitter::V_D_
private

Definition at line 74 of file LagrangeMultipliersFitter.h.

Referenced by applyLagrangianConstraints(), and computeVariance().

◆ weight_

double tauImpactParameter::LagrangeMultipliersFitter::weight_
private

Definition at line 66 of file LagrangeMultipliersFitter.h.

Referenced by setWeight().