Bug Summary

File:libraries/PID/DKinFit.cc
Location:line 933, column 15
Description:Value stored to 'numConstraints' is never read

Annotated Source Code

1#include <stdio.h>
2#include <stdlib.h>
3
4#include <cctype>
5#include <algorithm>
6#include "PID/DKinFit.h"
7using namespace std;
8/** @file DKinFit.C
9 * @brief DKinFit class source file.
10 */
11//_____________________________________________________________________________
12/** @class DKinFit
13 * @brief Utility class used to perform kinematic fitting.
14 *
15 * Specify the particles (i.e. DKinematicData) to be fit using SetInitial() and
16 * SetFinal() .
17 * Each DKinematicData must have a proper covariance matrix set.
18 *
19 * If needed you also SetMissingParticle() or SetMassConstraint()
20 * appropriately
21 *
22 * Then just Fit() !
23 *
24 * The results of the fit can be obtained thru various DKinFit member functions.
25 * All @a getter member functions return results of the last fit run.
26 *
27 */
28//_____________________________________________________________________________
29
30DKinFit::DKinFit(){
31 /// Default Constructor
32 _ndf = 0;
33 _chi2 = -1.;
34 _missingMass = -1.;
35 _missingParticle = false;
36 _extraC = 0;
37 _invariantMass = -1.;
38 _extraC_miss.clear();
39 _useTimingInformation = false;
40 _verbose = 0;
41}
42//_____________________________________________________________________________
43
44void DKinFit::_Copy(const DKinFit &__kfit){
45 /// Copy @a kfit data members.
46 _pulls = __kfit._pulls;
47 _chi2 = __kfit._chi2;
48 _ndf = __kfit._ndf;
49 _kDataInitial_in = __kfit._kDataInitial_in;
50 _kDataFinal_in = __kfit._kDataFinal_in;
51 _kDataInitial_out = __kfit._kDataInitial_out;
52 _kDataFinal_out = __kfit._kDataFinal_out;
53 _cov.ResizeTo(__kfit._cov);
54 _cov = __kfit._cov;
55 _missingMass = __kfit._missingMass;
56 _missingParticle = __kfit._missingParticle;
57 _extraC = __kfit._extraC;
58 _invariantMass = __kfit._invariantMass;
59 _extraC_meas = __kfit._extraC_meas;
60 _extraC_miss = __kfit._extraC_miss;
61 _useTimingInformation = __kfit._useTimingInformation;
62 _verbose = __kfit._verbose;
63 for(int i = 0; i < 3; i++) _sigma_missing[i] = __kfit._sigma_missing[i];
64}
65
66//_____________________________________________________________________________
67void DKinFit::Fit()
68{
69 /// Private function used internally to perform kinematic fitting.
70 _kDataInitial_out.clear();
71 _kDataFinal_out.clear();
72 _extraC_miss.clear();
73 const int numInitial = (int)_kDataInitial_in.size();
74 const int numFinal = (int)_kDataFinal_in.size();
75 const int numParts = numInitial + numFinal;
76 const int numypar = 3; // number of measured quantites
77 const int dim = numypar*numParts; // number of measured quantites
78 if(_verbose>0)
79 {
80 cerr << "numInitial: " << numInitial << endl;
81 cerr << "numFinal: " << numFinal << endl;
82 cerr << "numParts: " << numParts << endl;
83 cerr << "dim: " << dim << endl;
84 }
85 // Resize the covariance matrix
86 _cov.ResizeTo(dim, dim);
87
88 //const int dim = 3*numParts + 1; // number of measured quantites
89 // check to see that covariance matrix size is consistent with # particles
90 if(_cov.GetNrows() != dim || _cov.GetNcols() != dim)
91 {
92 // wrong size
93 std::cout << "Error! <DKinFit::_MainFitter> Covariance matrix size is NOT "
94 << "correct for current number of particles. For " << numParts
95 << " particles, the covariance matrix should be " << dim << " x "
96 << dim << ", but the matrix passed in is " << _cov.GetNrows()
97 << " x " << _cov.GetNcols() << std::endl;
98 abort();
99 }
100
101 for(int i=0;i<numInitial;i++)
102 {
103 _kDataInitial_out.push_back(_kDataInitial_in[i]); /// Make sure we're not using the same pointer
104 for(int j=0;j<numypar;j++)
105 {
106 for(int k=0;k<numypar;k++)
107 {
108 _cov(numypar*i+j,numypar*i+k) = (_kDataInitial_out[i].errorMatrix())(j,k);
109 }
110 }
111 }
112 for(int i=0;i<numFinal;i++)
113 {
114 // Convert from cartesian coordinates to the 5x1 state vector corresponding
115 // to the tracking error matrix
116 double vect[5];
117 DVector3 mom=_kDataFinal_in[i].momentum();
118 DVector3 pos=_kDataFinal_in[i].position();
119 double x=pos.x();
120 double y=pos.y();
121 vect[2]=tan(M_PI_21.57079632679489661923-mom.Theta());
122 vect[1]=mom.Phi();
123 double sinphi=sin(vect[1]);
124 double cosphi=cos(vect[1]);
125 vect[0]=_kDataFinal_in[i].charge()/mom.Perp();
126 vect[4]=pos.z();
127 vect[3]=pos.Perp();
128 if ((x>0 && sinphi>0) || (y<0 && cosphi>0) || (y>0 && cosphi<0)
129 || (x<0 && sinphi<0)) vect[3]*=-1.;
130
131 // Convert from the 5x5 tracking error matrix to the 7x7 matrix currently
132 // used by the kinematic fitter
133 _kDataFinal_in[i].setErrorMatrix(Get7x7ErrorMatrix(_kDataFinal_in[i].mass(),vect,_kDataFinal_in[i].TrackingErrorMatrix()));
134
135
136 _kDataFinal_out.push_back(_kDataFinal_in[i]); /// Make sure we're not using the same pointer
137 for(int j=0;j<numypar;j++)
138 {
139 for(int k=0;k<numypar;k++)
140 {
141 _cov( numypar*i+j + numypar*numInitial, numypar*i+k + numypar*numInitial) = (_kDataFinal_out[i].errorMatrix())(j,k);
142 }
143 }
144 }
145
146 if(_verbose>1)
147 {
148 cerr << "Total error matrix: " << endl;
149 for(int j=0;j<_cov.GetNrows();j++)
150 {
151 for(int k=0;k<_cov.GetNcols();k++)
152 {
153 cerr << _cov(j,k) << " ";
154 }
155 cerr << endl;
156 }
157 }
158 // Chi2 of fit
159 double chi2=1e16,chi2old=0.;
160
161
162 // Get the number of constraint equations
163 int numConstraints = 4;
164 if(_extraC) numConstraints += _extraC;
165 // For Two gammas
166 //numConstraints = 4;
167
168 // get the degrees of freedom
169 _ndf = 4;
170 if(_missingParticle) _ndf -= 3;
171 _ndf += _extraC;
172 //if(_extraC) _ndf += 1;
173 // For Two gammas
174 //_ndf = 1;
175
176 if(_verbose>0) cerr << "numConstraints/ndf: " << numConstraints << " " << _ndf << endl;
177
178 int i;
179 // The following were changed to vectors to avoid compiler warnings on SunOS 11/12/2007 D.L.
180 vector<double> mass(numParts); //double mass[numParts];
181 vector<double> energy(numParts); //double energy[numParts];
182 double missingEnergy = 999.0;
183 vector<int> initialOrFinal(numParts); //int initialOrFinal[numParts];
184 // The following 3 lines were commented out to avoid compiler warnings 7/31/07 D.L.
185 //double p[numParts],erg[numParts],px[numParts],py[numParts],pz[numParts];
186 //double e_miss=0.,e_inv=0.,px_inv=0.,py_inv=0.,pz_inv=0.;
187 //double dumt, dumx, dumy, dumz;
188 DVector3 p3_miss;
189 // The following were changed to vectors to avoid compiler warnings on SunOS 11/12/2007 D.L.
190 vector<double> pxMeas(_extraC), pyMeas(_extraC), pzMeas(_extraC), EMeas(_extraC);
191 vector<double> pxTot(_extraC), pyTot(_extraC), pzTot(_extraC), ETot(_extraC), massTot(_extraC);
192 //double pxMeas[_extraC], pyMeas[_extraC], pzMeas[_extraC], EMeas[_extraC];
193 //double pxTot[_extraC], pyTot[_extraC], pzTot[_extraC], ETot[_extraC], massTot[_extraC];
194
195 for(int i=0;i<numParts;i++)
196 {
197 if(i<numInitial) initialOrFinal[i] = +1;
198 else initialOrFinal[i] = -1;
199 }
200
201 /*********** Define matricies needed to perform kinematic fitting **********/
202
203 // Note: for any matrix m, mT is the transpose of m.
204 //
205 DLorentzVector missingMomentum;
206
207 // y(i) ==> measured particle quantities.
208 // i = 3*n+0,3*n+1,3*n+2 are particle n's p_x,p_y,p_z respectively.
209 DMatrix yi(dim,1),y(dim,1); // yi will be intial values
210 // x(i) ==> missing particle quantities (0,1,2) = (px,py,pz)
211 DMatrix x(3,1);
212 //DMatrix x(4,1);
213 // a(i,j) = dc(i)/dx(j) ==> derivatives of constraint eq's w/r to missing
214 // particle quantities x.
215 DMatrix a(numConstraints,3),aT(3,numConstraints);
216 //DMatrix a(numConstraints,4),aT(4,numConstraints);
217 // b(i,j) = dc(i)/dy(j) ==> derivatives of constraint eq's w/r to measured
218 // quantities y.
219 DMatrix b(numConstraints,dim),bT(dim,numConstraints);
220 // c(i) is constraint eq. i ==> (0,1,2,3) are (e,px,py,pz) constraints, any
221 // extra mass constraints are i = 4 and up.
222 DMatrix c(numConstraints,1);
223 DMatrix c0(numConstraints,1);
224 // eps(i) = yi(i)-y(i) are the overall changes in the measured quantities
225 DMatrix eps(dim,1),epsT(1,dim);
226 // delta(i) is used to update the measured quantities. After an iteration of
227 // the fitting process, y =yi-delta is how y is updated.
228 DMatrix delta(dim,1);
229 // xsi(i) same as delta (above) but for x.
230 DMatrix xsi(3,1);
231 //DMatrix xsi(4,1);
232 // gb is a utility matrix used
233 DMatrix gb(numConstraints,numConstraints);
234 DMatrix gb0(numConstraints,numConstraints);
235 // cov_fit is the covariance matrix OUTPUT by the fit.
236 DMatrix cov_fit(dim,dim);
237 // Lagrange multiplier array
238 DMatrix L(numConstraints,1);
239 DMatrix LT(1,numConstraints);
240
241 /***************************** Initial Set Up ******************************/
242
243 // set up vector of initial measured values:
244 for(i = 0; i < numInitial; i++)
245 {
246 mass[i] = _kDataInitial_in[i].mass(); // particle's mass
247 yi(numypar*i + 0,0) = _kDataInitial_in[i].px();
248 yi(numypar*i + 1,0) = _kDataInitial_in[i].py();
249 yi(numypar*i + 2,0) = _kDataInitial_in[i].pz();
250 if(_verbose) cerr << "initial masses: " << mass[i] << endl;
251 }
252 for(i = 0; i < numFinal; i++)
253 {
254 mass[numInitial + i] = _kDataFinal_in[i].mass(); // particle's mass
255 yi(numypar*(numInitial+i) + 0,0) = _kDataFinal_in[i].px();
256 yi(numypar*(numInitial+i) + 1,0) = _kDataFinal_in[i].py();
257 yi(numypar*(numInitial+i) + 2,0) = _kDataFinal_in[i].pz();
258 if(_verbose) cerr << "final masses: " << mass[i+numInitial] << endl;
259 }
260 y = yi; // start off y to be yi
261
262 // For Two gammas
263 if(_missingParticle)
264 {
265 x(0,0) = 0.0;
266 x(1,0) = 0.0;
267 x(2,0) = 0.0;
268 for(i = 0; i < numInitial; i++)
269 {
270 x(0,0) += _kDataInitial_in[i].px();
271 x(1,0) += _kDataInitial_in[i].py();
272 x(2,0) += _kDataInitial_in[i].pz();
273 }
274 for(i = 0; i < numFinal; i++)
275 {
276 x(0,0) -= _kDataFinal_in[i].px();
277 x(1,0) -= _kDataFinal_in[i].py();
278 x(2,0) -= _kDataFinal_in[i].pz();
279 }
280 }
281
282 /************************* Begin Iteration Process *************************/
283 for(int iter = 0; iter < 10; iter++)
284 {
285 // it should converge by 10 iterations
286 // reset all utility matricies
287 a.Zero();
288 aT.Zero();
289 b.Zero();
290 bT.Zero();
291 c.Zero();
292 delta.Zero();
293 xsi.Zero();
294
295 if(_verbose>1)
296 {
297 cerr << "y: " << endl;
298 for(int j=0;j<y.GetNrows();j++)
299 {
300 for(int k=0;k<y.GetNcols();k++)
301 {
302 cerr << j << " " << k << " " << y(j,k) << " ";
303 }
304 cerr << endl;
305 }
306 }
307
308 if(_verbose>1)
309 {
310 cerr << "x: " << endl;
311 for(int j=0;j<x.GetNrows();j++)
312 {
313 for(int k=0;k<x.GetNcols();k++)
314 {
315 cerr << j << " " << k << " " << x(j,k) << " ";
316 }
317 cerr << endl;
318 }
319 }
320
321 // Set up the constraint eq's for this iteration:
322 // c(0) (energy constraint): e_f - e_i = 0
323 // c(i) (p3 constraint): p3_f - p3_i = 0
324 if(_verbose>1) cerr << "setting up the constraints " << endl;
325 for(int j=0;j<4;j++) c0(j,0) = 0.;
326 for(int k=0;k<numInitial;k++)
327 {
328 int offset = k*3;
329 energy[k] = sqrt(pow(mass[k],2) + pow(y(offset+0,0),2) + pow(y(offset+1,0),2) + pow(y(offset+2,0),2));
330 c0(0,0) += y(offset+0,0);
331 c0(1,0) += y(offset+1,0);
332 c0(2,0) += y(offset+2,0);
333 c0(3,0) += energy[k];
334 }
335 for(int k=0;k<numFinal;k++)
336 {
337 int offset = (k+numInitial)*3;
338 energy[k+numInitial] = sqrt(pow(mass[k+numInitial],2) + pow(y(offset+0,0),2) + pow(y(offset+1,0),2) + pow(y(offset+2,0),2));
339 c0(0,0) -= y(offset+0,0);
340 c0(1,0) -= y(offset+1,0);
341 c0(2,0) -= y(offset+2,0);
342 c0(3,0) -= energy[k+numInitial];
343 }
344
345 if(_missingParticle)
346 {
347 missingEnergy = sqrt(pow(_missingMass,2) + pow(x(0,0),2) + pow(x(1,0),2) + pow(x(2,0),2));
348 c0(0,0) -= x(0,0);
349 c0(1,0) -= x(1,0);
350 c0(2,0) -= x(2,0);
351 c0(3,0) -= missingEnergy;
352 }
353
354 /// If there is a extra mass constraints
355 for(int j=0;j<_extraC;j++)
356 {
357 int cindex = 4 + j;
358 int yoffset = 3*numInitial;
359 int numConstraintParticles = (int)_constraintParticles[j].size();
360 pxMeas[j]=pyMeas[j]=pzMeas[j]=EMeas[j]=0.0;
361 pxTot[j]=pyTot[j]=pzTot[j]=ETot[j]=0.0;
362 massTot[j] = 0.0;
363 _extraC_miss.push_back(false); // Assume missing particle is not in constraint
364 /// Add in the particles to the constraint
365 for(int k=0;k<numConstraintParticles;k++)
366 {
367 int partIndex = _constraintParticles[j][k];
368 int Eindex = numInitial + partIndex;
369 /// Use a measured particle
370 if(partIndex>=0 && partIndex<numFinal)
371 {
372 if(_verbose>1) cerr << "Using constraint final particle: " << partIndex << endl;
373 pxMeas[j] += y(yoffset + 3*partIndex + 0 ,0);
374 pyMeas[j] += y(yoffset + 3*partIndex + 1 ,0);
375 pzMeas[j] += y(yoffset + 3*partIndex + 2 ,0);
376 EMeas[j] += energy[Eindex];
377 if(_verbose>1) cerr << "constraint4vec: " << pxMeas[j] << " " << pyMeas[j] << " " << pzMeas[j] << " " << EMeas[j] << endl;
378 }
379 else if(partIndex == -1)
380 {
381 _extraC_miss[j] = true; // Assume missing particle is not in constraint
382 }
383 }
384 pxTot[j] = pxMeas[j];
385 pyTot[j] = pyMeas[j];
386 pzTot[j] = pzMeas[j];
387 ETot[j] = EMeas[j];
388 if(_extraC_miss[j])
389 {
390 pxTot[j] += x(0,0);
391 pyTot[j] += x(1,0);
392 pzTot[j] += x(2,0);
393 ETot[j] += missingEnergy;
394 }
395 massTot[j] = sqrt(pow(ETot[j],2) - pow(pxTot[j],2) - pow(pyTot[j],2) - pow(pzTot[j],2));
396 if(_verbose>1) cerr << "massTot/constraintMass: " << massTot[j] << " " << _constraintMasses[j] << endl;
397 c0(cindex, 0) = massTot[j] - _constraintMasses[j];
398 }
399
400 if(_verbose)
401 {
402 for(int k=0;k<numParts;k++) cerr << "energy: " << energy[k] << endl;
403 cerr << "missing energy: " << missingEnergy << endl;
404 }
405
406 // This is for when I have 4 constraints
407 //float gammaE0 = sqrt(y(0,0)*y(0,0) + y(1,0)*y(1,0) + y(2,0)*y(2,0));
408 //float gammaE1 = sqrt(y(3,0)*y(3,0) + y(4,0)*y(4,0) + y(5,0)*y(5,0));
409 //float pi0E = sqrt(x(0,0)*x(0,0) + x(1,0)*x(1,0) + x(2,0)*x(2,0) + _missingMass * _missingMass);
410 //c(0,0) = x(0,0) - y(0,0) - y(3,0);
411 //c(1,0) = x(1,0) - y(1,0) - y(4,0);
412 //c(2,0) = x(2,0) - y(2,0) - y(5,0);
413 //c(3,0) = pi0E - gammaE0 - gammaE1;
414 //if(_verbose>1) cerr << "set up the constraints " << endl;
415
416 if(_verbose>1)
417 {
418 cerr << "c: " << endl;
419 for(int j=0;j<c0.GetNrows();j++)
420 {
421 for(int k=0;k<c0.GetNcols();k++)
422 {
423 cerr << j << " " << k << " " << c0(j,k) << " ";
424 }
425 cerr << endl;
426 }
427 }
428
429
430 //This is for 4 constraints
431 // Derivatives of constraint equations wrt unknown quantities.
432 if(_missingParticle) /// NEED TO DO THIS RIGHT!!!!!!!!!!!!! MB
433 {
434 for(int j=0;j<3;j++)
435 {
436 a(j,j) = -1.0; // derivatives of p3 constraint eq's
437 a(3,j) = -x(j,0)/missingEnergy; // derivatives of energy constraint eq.
438 }
439 for(int j=0;j<_extraC;j++)
440 {
441 int aindex = 4 + j;
442 if(_extraC_miss[j])
443 {
444 a(aindex, 0) = (((x(0,0)*ETot[j])/missingEnergy) - pxTot[j])/massTot[j];
445 a(aindex, 1) = (((x(1,0)*ETot[j])/missingEnergy) - pyTot[j])/massTot[j];
446 a(aindex, 2) = (((x(2,0)*ETot[j])/missingEnergy) - pzTot[j])/massTot[j];
447 }
448 }
449 }
450
451 if(_verbose>1)
452 {
453 cerr << "a: " << endl;
454 for(int j=0;j<a.GetNrows();j++)
455 {
456 for(int k=0;k<a.GetNcols();k++)
457 {
458 cerr << a(j,k) << " ";
459 }
460 cerr << endl;
461 }
462 }
463
464 aT.Transpose(a); // set a transpose
465
466 if(_verbose>1)
467 {
468 cerr << "aT: " << endl;
469 for(int j=0;j<aT.GetNrows();j++)
470 {
471 for(int k=0;k<aT.GetNcols();k++)
472 {
473 cerr << aT(j,k) << " ";
474 }
475 cerr << endl;
476 }
477 }
478
479 // Set up derivatives matrix b(i,j) = dc(i)/dy(j):
480 //This is for 4 constraints
481 for(int k=0;k<numInitial;k++)
482 {
483 int offset = k*3;
484 b(0,offset+0) = 1.0;
485 b(1,offset+1) = 1.0;
486 b(2,offset+2) = 1.0;
487 for(int j=0;j<3;j++) b(3,offset+j) = y(offset+j,0)/energy[k];
488 }
489 for(int k=0;k<numFinal;k++)
490 {
491 int offset = (k+numInitial)*3;
492 b(0,offset+0) = -1.0;
493 b(1,offset+1) = -1.0;
494 b(2,offset+2) = -1.0;
495 for(int j=0;j<3;j++) b(3,offset+j) = -y(offset+j,0)/energy[k+numInitial];
496 }
497 for(int j=0;j<_extraC;j++)
498 {
499 int bindex = 4 + j;
500 int yoffset = 3*numInitial;
501 int numConstraintParticles = (int)_constraintParticles[j].size();
502 /// Add in the particles to the constraint
503 for(int k=0;k<numConstraintParticles;k++)
504 {
505 int partIndex = _constraintParticles[j][k];
506 int Eindex = numInitial + partIndex;
507 /// Use a measured particle
508 if(partIndex>=0 && partIndex<numFinal)
509 {
510 b(bindex, yoffset + 3*partIndex +0) = (((y(yoffset + 3*partIndex +0,0)*ETot[j])/energy[Eindex]) - pxTot[j])/massTot[j];
511 b(bindex, yoffset + 3*partIndex +1) = (((y(yoffset + 3*partIndex +1,0)*ETot[j])/energy[Eindex]) - pyTot[j])/massTot[j];
512 b(bindex, yoffset + 3*partIndex +2) = (((y(yoffset + 3*partIndex +2,0)*ETot[j])/energy[Eindex]) - pzTot[j])/massTot[j];
513 }
514 }
515 }
516
517 c=c0+b*(yi-y);
518
519
520 if(_verbose>1)
521 {
522 cerr << "b: " << endl;
523 for(int j=0;j<b.GetNrows();j++)
524 {
525 for(int k=0;k<b.GetNcols();k++)
526 {
527 cerr << b(j,k) << " ";
528 }
529 cerr << endl;
530 }
531 }
532
533 bT.Transpose(b); // set b transpose
534
535 if(_verbose>1)
536 {
537 cerr << "bT: " << endl;
538 for(int j=0;j<bT.GetNrows();j++)
539 {
540 for(int k=0;k<bT.GetNcols();k++)
541 {
542 cerr << bT(j,k) << " ";
543 }
544 cerr << endl;
545 }
546 }
547
548 // gb is a utility matrix we'll need a lot below
549 gb = b * _cov * bT;
550 gb0=gb;
551
552 if(_verbose>1)
553 {
554 cerr << "gb: " << endl;
555 for(int j=0;j<gb.GetNrows();j++)
556 {
557 for(int k=0;k<gb.GetNcols();k++)
558 {
559 cerr << gb(j,k) << " ";
560 }
561 cerr << endl;
562 }
563 }
564
565 if( gb.Determinant() == 0 ) break; // for the case when gb is singular
566
567 gb.Invert();
568 if(_verbose>1)
569 {
570 cerr << "gb inverted: " << endl;
571 for(int j=0;j<gb.GetNrows();j++)
572 {
573 for(int k=0;k<gb.GetNcols();k++)
574 {
575 cerr << gb(j,k) << " ";
576 }
577 cerr << endl;
578 }
579
580 }
581
582 if(_missingParticle)
583 {
584 if(_verbose>1)
585 {
586 cerr << "(aT * gb * a): " << endl;
587 for(int j=0;j<(aT * gb * a).GetNrows();j++)
588 {
589 for(int k=0;k<(aT * gb * a).GetNcols();k++)
590 {
591 cerr << (aT * gb * a)(j,k) << " ";
592 }
593 cerr << endl;
594 }
595 cerr << "Before invert: " << endl;
596 cerr << "Determinant: " << (aT * gb * a).Determinant() << endl;
597 }
598 //(aT * gb * a).Invert();
599 if(_verbose>1) cerr << "After invert: " << endl;
600 // Is this what we should do?
601 if( (aT * gb * a).Determinant() == 0.0 ) break;
602
603 // calculate the change to be made to the missing particle quantities
604 xsi = ((aT * gb * a).Invert())*(aT * gb * c);
605 // update x
606 x -= xsi;
607 // calculate the changes to be made to the measured quantities
608 delta = _cov * bT * gb * (c - a*xsi);
609
610 if(_verbose>1)
611 {
612 cerr << "xsi: " << endl;
613 for(int j=0;j<xsi.GetNrows();j++)
614 {
615 for(int k=0;k<xsi.GetNcols();k++)
616 {
617 cerr << xsi(j,k) << " ";
618 }
619 cerr << endl;
620 }
621
622 cerr << "(c - a*xsi): " << endl;
623 for(int j=0;j<(c - a*xsi).GetNrows();j++)
624 {
625 for(int k=0;k<(c - a*xsi).GetNcols();k++)
626 {
627 cerr << (c - a*xsi)(j,k) << " ";
628 }
629 cerr << endl;
630 }
631
632 }
633
634 }
635 else
636 {
637 delta = _cov * bT * gb * c;
638 }
639
640 // update the measured quantities y
641 //y -= delta;
642 y=yi-delta;
643
644 if(_verbose>1)
645 {
646 cerr << "delta: " << endl;
647 for(int j=0;j<delta.GetNrows();j++)
648 {
649 for(int k=0;k<delta.GetNcols();k++)
650 {
651 cerr << delta(j,k) << " ";
652 }
653 cerr << endl;
654 }
655 }
656
657 if(_verbose>1)
658 {
659 cerr << "new y: " << endl;
660 for(int j=0;j<y.GetNrows();j++)
661 {
662 for(int k=0;k<y.GetNcols();k++)
663 {
664 cerr << y(j,k) << " ";
665 }
666 cerr << endl;
667 }
668 }
669
670
671 eps = yi - y;
672 epsT.Transpose(eps);
673
674 // Lagrange multiplier array and its transpose
675 L=gb*c;
676 LT.Transpose(L);
677
678 // Compute current chi2 value and check for convergence
679 chi2=(LT*gb0*L+2.*LT*c0)(0,0);
680 if (fabs(chi2-chi2old)<0.001) break;
681 chi2old=chi2;
682
683 if(_verbose>1)
684 {
685 cerr << "eps: " << endl;
686 for(int j=0;j<eps.GetNrows();j++)
687 {
688 for(int k=0;k<eps.GetNcols();k++)
689 {
690 cerr << eps(j,k) << " ";
691 }
692 cerr << endl;
693 }
694 }
695
696
697 } /* end iteration process */
698
699 // get the total changes in the measured quantities
700 eps = yi - y;
701 epsT.Transpose(eps);
702
703 // get the fit covariance matrix
704 if(_missingParticle)
705 {
706 if( (aT * gb * a).Determinant() != 0.0 ) {
707 cov_fit = _cov - _cov*bT*gb*b*_cov + _cov*bT*gb*a*((aT*gb*a).Invert())*aT*gb*b*_cov;
708 }
709 else
710 {
711 if(_verbose>1) cerr << "Not going to fill the missing cov matrix output because of singularity...." << endl;
712 }
713 }
714 else
715 {
716 cov_fit = _cov - _cov*bT*gb*b*_cov;
717 }
718
719 DMatrixDSym dumerr(7);
720 //dumerr(0,0) = 1.0;
721 //cerr << "FILLING ERROR MATRIX" << endl;
722 //cerr << "dum row/col: " << dumerr.GetNrows() << " " << dumerr.GetNcols() << endl;
723 //cerr << "kin row/col: " << _kDataInitial_out[0].errorMatrix().GetNrows() << " " << _kDataInitial_out[0].errorMatrix().GetNcols() << endl;
724 //_kDataInitial_out[0].setErrorMatrix(dumerr);
725 //cerr << "FILLED ERROR MATRIX" << endl;
726
727 for(int i=0;i<numInitial;i++)
728 {
729 dumerr = _kDataInitial_out[i].errorMatrix();
730 for(int j=0;j<numypar;j++)
731 {
732 for(int k=0;k<numypar;k++)
733 {
734 dumerr(j,k) = cov_fit(numypar*i+j,numypar*i+k);
735 }
736 }
737 _kDataInitial_out[i].setErrorMatrix(dumerr);
738 }
739 for(int i=0;i<numFinal;i++)
740 {
741 dumerr = _kDataFinal_out[i].errorMatrix();
742 for(int j=0;j<numypar;j++)
743 {
744 for(int k=0;k<numypar;k++)
745 {
746 dumerr(j,k) = cov_fit( numypar*i+j + numypar*numInitial, numypar*i+k + numypar*numInitial);
747 }
748 }
749 _kDataFinal_out[i].setErrorMatrix(dumerr);
750 }
751
752 // get chi2 and the pulls
753 _pulls.resize(dim);
754 for(i = 0; i < dim; i++)
755 {
756 _pulls[i] = -eps(i,0)/sqrt(_cov(i,i) - cov_fit(i,i));
757 }
758
759 // _chi2 = (epsT * bT * gb * b * eps)(0,0); // the (0,0)...only...element
760 _chi2=chi2;
761
762
763 // DRandom rnd;
764 // rnd.SetSeed((unsigned int)(10000*_chi2));
765
766
767 if(_verbose>0) cerr << "Prob: " << TMath::Prob(_chi2, _ndf) << " " << _chi2 << " " << _ndf << endl;
768
769 // update output of measured quantities
770 if(_verbose>1) cerr << "Filling the 4vecs output...." << endl;
771 for(i = 0; i < numInitial; i++)
772 {
773 int offset = 3*i;
774 _kDataInitial_out[i].setMass(mass[i]);
775 _kDataInitial_out[i].setMomentum(DVector3(y(offset+0,0), y(offset+1,0), y(offset+2,0)));
776 }
777 for(i = 0; i < numFinal; i++)
778 {
779 int offset = 3*(i+numInitial);
780 _kDataFinal_out[i].setMass(mass[i+numInitial]);
781 _kDataFinal_out[i].setMomentum(DVector3(y(offset+0,0), y(offset+1,0), y(offset+2,0)));
782 }
783 if(_verbose>1) cerr << "Filled the 4vecs output...." << endl;
784
785 if(_useTimingInformation)
786 {
787 //_chi2 = 0;
788 //_ndf = 0;
789 for(i = 0; i < numFinal; i++)
790 {
791 if(_verbose>=1) cerr << "detector: " << i << " " << _kDataFinal_in[i].t1_detector() << endl;
792 //if(_kDataFinal_in[i].t1_detector() != SYS_NULL)
793 //if(_kDataFinal_in[i].t1_detector() == SYS_TOF)
794 //if(_kDataFinal_in[i].t1_detector() == SYS_BCAL)
795 if(_kDataFinal_in[i].t1_detector() == SYS_BCAL || _kDataFinal_in[i].t1_detector() == SYS_TOF)
796 {
797 //cerr << "Before: chi2/_ndf/prob: " << _chi2 << " " << _ndf << " " << TMath::Prob(_chi2, _ndf) << endl;
798 double beta = _kDataFinal_out[i].lorentzMomentum().Beta();
799 double deltaBeta = _kDataFinal_out[i].deltaBeta();
800 double measuredBeta_err = _kDataFinal_out[i].measuredBeta_err();
801 double px = _kDataFinal_out[i].lorentzMomentum().X();
802 double py = _kDataFinal_out[i].lorentzMomentum().Y();
803 double pz = _kDataFinal_out[i].lorentzMomentum().Z();
804 double pxerr = sqrt(_kDataFinal_out[i].errorMatrix()(0,0));
805 double pyerr = sqrt(_kDataFinal_out[i].errorMatrix()(1,1));
806 double pzerr = sqrt(_kDataFinal_out[i].errorMatrix()(2,2));
807 double X = beta*beta;
808 double A = pow(_kDataFinal_out[i].lorentzMomentum().Rho(),2);
809 double B = pow(_kDataFinal_out[i].lorentzMomentum().E(),2);
810 double errA = sqrt(pow(px*2.0*pxerr,2) + pow(py*2.0*pyerr,2) + pow(pz*2.0*pzerr,2) );
811 double errB = errA;
812 double errX = X*sqrt( pow(errA/A,2) + pow(errB/B,2) - 2.0*errA*errB/(A*B));
813 double calcBeta_err = beta*(1.0/2.0)*( errX/X );
814 double timing_chi2 = pow(deltaBeta,2)/(pow(calcBeta_err,2) + pow(measuredBeta_err,2));
815 //double timing_chi2 = pow(deltaBeta,2)/(pow(calcBeta_err,2));
816 if(_verbose>=1) cerr << "deltaBeta/err/timing_chi2: " << i << " " << deltaBeta << " " << measuredBeta_err << " " << timing_chi2 << endl;
817 if(_verbose>=1) cerr << "Beta/measured: " << i << " " << _kDataFinal_out[i].lorentzMomentum().Beta() << "\t";
818 if(_verbose>=1) cerr << _kDataFinal_out[i].measuredBeta() << endl;
819 if(_verbose>=1) cerr << "t1_detector: " << i << " " << _kDataFinal_in[i].t1_detector() << "\t";
820 if(_verbose>=1) cerr << "chi2: " << _chi2 << endl;
821 _chi2 += timing_chi2;
822 _ndf += 1;
823 //cerr << "chi2/_ndf/prob: " << _chi2 << " " << _ndf << " " << TMath::Prob(_chi2, _ndf) << endl;
824 }
825 }
826 }
827 // set the missing particle covariance matrix
828 if(_verbose>1) cerr << "Filling the missing cov matrix output...." << endl;
829 // Not sure about this part
830 if(_missingParticle)
831 {
832 if( (aT * gb * a).Determinant() != 0.0 )
833 {
834 if(_missingParticle)
835 this->_SetMissingParticleErrors((aT * gb * a).Invert(),x);
836 if(_verbose>1) cerr << "Filled the missing cov matrix output...." << endl;
837 }
838 else
839 {
840 if(_verbose>1) cerr << "Not going to fill the missing cov matrix output 'cos singular...." << endl;
841 }
842 }
843}
844//_____________________________________________________________________________
845
846void DKinFit::FitTwoGammas(const float __missingMass, const float errmatrixweight=1.0)
847{
848 /// Private function used internally to perform kinematic fitting.
849 //const int numParts = (int)_p4in.size();
850 //int dum_dim = 0;
851 //for(int i=0;i<(int)_kDataInitial_in.size();i++)
852 //{
853 //dum_numParts += _kDataInitial_in[i].NumQuantitiesForDKinFit();
854 //}
855 _missingMass = __missingMass;
856 _kDataInitial_out.clear();
857 _kDataFinal_out.clear();
858 //float errmatrixweight = 1.0;
859 const int numInitial = (int)_kDataInitial_in.size();
860 const int numFinal = (int)_kDataFinal_in.size();
861 const int numParts = numInitial + numFinal;
862 //const int dim = 4*numParts; // number of measured quantites
863 const int numypar = 3; // number of measured quantites
864 const int dim = numypar*numParts; // number of measured quantites
865 if(_verbose>0)
866 {
867 cerr << "numInitial: " << numInitial << endl;
868 cerr << "numFinal: " << numFinal << endl;
869 cerr << "numParts: " << numParts << endl;
870 cerr << "dim: " << dim << endl;
871 }
872 // Resize the covariance matrix
873 _cov.ResizeTo(dim, dim);
874
875 //const int dim = 3*numParts + 1; // number of measured quantites
876 // check to see that covariance matrix size is consistent with # particles
877 if(_cov.GetNrows() != dim || _cov.GetNcols() != dim)
878 {
879 // wrong size
880 std::cout << "Error! <DKinFit::_MainFitter> Covariance matrix size is NOT "
881 << "correct for current number of particles. For " << numParts
882 << " particles, the covariance matrix should be " << dim << " x "
883 << dim << ", but the matrix passed in is " << _cov.GetNrows()
884 << " x " << _cov.GetNcols() << std::endl;
885 abort();
886 }
887
888 //if(_verbose>1)
889 {
890 for(int i=0;i<(int)_kDataInitial_in.size();i++)
891 {
892 _kDataInitial_out.push_back(_kDataInitial_in[i]); /// Make sure we're not using the same pointer
893 for(int j=0;j<numypar;j++)
894 {
895 for(int k=0;k<numypar;k++)
896 {
897 _cov(numypar*i+j,numypar*i+k) = (_kDataInitial_out[i].errorMatrix())(j,k)/errmatrixweight;
898 }
899 }
900 }
901 for(int i=0;i<(int)_kDataFinal_in.size();i++)
902 {
903 _kDataFinal_out.push_back(_kDataFinal_in[i]); /// Make sure we're not using the same pointer
904 for(int j=0;j<numypar;j++)
905 {
906 for(int k=0;k<numypar;k++)
907 {
908 _cov(numypar*i+j,numypar*i+k) = (_kDataFinal_out[i].errorMatrix())(j,k)/errmatrixweight;
909 }
910 }
911 }
912 }
913
914 if(_verbose>1)
915 {
916 cerr << "Total error matrix: " << endl;
917 for(int j=0;j<_cov.GetNrows();j++)
918 {
919 for(int k=0;k<_cov.GetNcols();k++)
920 {
921 cerr << _cov(j,k) << " ";
922 }
923 cerr << endl;
924 }
925 }
926
927
928 // The following was commented out to avoid compiler warnings 7/31/07 D.L.
929 //bool isEnergyMeasured[numInitial + numFinal];
930
931 // Get the number of constraint equations
932 int numConstraints = 4;
933 if(_extraC) numConstraints = 5;
Value stored to 'numConstraints' is never read
934 // For Two gammas
935 numConstraints = 4;
936
937 // get the degrees of freedom
938 _ndf = 4;
939 if(_missingParticle) _ndf -= 3;
940 if(_extraC) _ndf += 1;
941 // For Two gammas
942 _ndf = 1;
943
944 if(_verbose>0) cerr << "numConstraints/ndf: " << numConstraints << " " << _ndf << endl;
945
946 int i;
947 // The following were changed to vectors to avoid compiler warnings on SunOS 11/12/2007 D.L.
948 vector<double> mass(numParts); //double mass[numParts];
949 vector<int> initialOrFinal(numParts); //int initialOrFinal[numParts];
950
951 // The following 3 lines were commented out to avoid compiler warnings 7/31/07 D.L.
952 //double p[numParts],erg[numParts],px[numParts],py[numParts],pz[numParts];
953 //double e_miss=0.,e_inv=0.,px_inv=0.,py_inv=0.,pz_inv=0.;
954 //double dumt, dumx, dumy, dumz;
955 DVector3 p3_miss;
956
957 for(int i=0;i<numParts;i++)
958 {
959 if(i<numInitial) initialOrFinal[i] = +1;
960 else initialOrFinal[i] = -1;
961 }
962
963 /*********** Define matricies needed to perform kinematic fitting **********/
964
965 // Note: for any matrix m, mT is the transpose of m.
966 //
967 DLorentzVector missingMomentum;
968
969 // y(i) ==> measured particle quantities.
970 // i = 3*n+0,3*n+1,3*n+=2 are particle n's p_x,p_y,p_z respectively.
971 DMatrix yi(dim,1),y(dim,1); // yi will be intial values
972 // x(i) ==> missing particle quantities (0,1,2) = (px,py,pz)
973 DMatrix x(3,1);
974 //DMatrix x(4,1);
975 // a(i,j) = dc(i)/dx(j) ==> derivatives of constraint eq's w/r to missing
976 // particle quantities x.
977 DMatrix a(numConstraints,3),aT(3,numConstraints);
978 //DMatrix a(numConstraints,4),aT(4,numConstraints);
979 // b(i,j) = dc(i)/dy(j) ==> derivatives of constraint eq's w/r to measured
980 // quantities y.
981 DMatrix b(numConstraints,dim),bT(dim,numConstraints);
982 // c(i) is constraint eq. i ==> (0,1,2,3) are (e,px,py,pz) constraints, any
983 // extra mass constraints are i = 4 and up.
984 DMatrix c(numConstraints,1);
985 // eps(i) = yi(i)-y(i) are the overall changes in the measured quantities
986 DMatrix eps(dim,1),epsT(1,dim);
987 // delta(i) is used to update the measured quantities. After an iteration of
988 // the fitting process, y -= delta is how y is updated.
989 DMatrix delta(dim,1);
990 // xsi(i) same as delta (above) but for x.
991 DMatrix xsi(3,1);
992 //DMatrix xsi(4,1);
993 // gb is a utility matrix used
994 DMatrix gb(numConstraints,numConstraints);
995 // cov_fit is the covariance matrix OUTPUT by the fit.
996 DMatrix cov_fit(dim,dim);
997
998 /***************************** Initial Set Up ******************************/
999
1000 // set up vector of initial measured values:
1001 for(i = 0; i < numInitial; i++)
1002 {
1003 mass[i] = _kDataInitial_in[i].mass(); // particle's mass
1004 yi(numypar*i + 0,0) = _kDataInitial_in[i].px();
1005 yi(numypar*i + 1,0) = _kDataInitial_in[i].py();
1006 yi(numypar*i + 2,0) = _kDataInitial_in[i].pz();
1007 }
1008 for(i = 0; i < numFinal; i++)
1009 {
1010 mass[numInitial + i] = _kDataFinal_in[i].mass(); // particle's mass
1011 yi(numypar*(numInitial+i) + 0,0) = _kDataFinal_in[i].px();
1012 yi(numypar*(numInitial+i) + 1,0) = _kDataFinal_in[i].py();
1013 yi(numypar*(numInitial+i) + 2,0) = _kDataFinal_in[i].pz();
1014 }
1015 y = yi; // start off y to be yi
1016
1017 // For Two gammas
1018 //if(_missingParticle)
1019 {
1020 //missingMomentum = this->MissingMomentum(_kDataInitial_in, _kDataFinal_in);
1021 x(0,0) = yi(0,0) + yi(3,0);
1022 x(1,0) = yi(1,0) + yi(4,0);
1023 x(2,0) = yi(2,0) + yi(5,0);
1024 }
1025
1026 /************************* Begin Iteration Process *************************/
1027 for(int iter = 0; iter < 10; iter++)
1028 {
1029 // it should converge by 10 iterations
1030 // reset all utility matricies
1031 a.Zero();
1032 aT.Zero();
1033 b.Zero();
1034 bT.Zero();
1035 c.Zero();
1036 delta.Zero();
1037 xsi.Zero();
1038
1039 // Set up the constraint eq's for this iteration:
1040 // c(0) (energy constraint): e_f - e_i = 0
1041 // c(i) (p3 constraint): p3_f - p3_i = 0
1042 for(int j=0;j<numConstraints;j++) c(j,0) = 0.;
1043
1044 // This is for when I have 4 constraints
1045 if(_verbose>1) cerr << "setting up the constraints " << endl;
1046 float gammaE0 = sqrt(y(0,0)*y(0,0) + y(1,0)*y(1,0) + y(2,0)*y(2,0));
1047 float gammaE1 = sqrt(y(3,0)*y(3,0) + y(4,0)*y(4,0) + y(5,0)*y(5,0));
1048 float pi0E = sqrt(x(0,0)*x(0,0) + x(1,0)*x(1,0) + x(2,0)*x(2,0) + __missingMass * __missingMass);
1049 c(0,0) = x(0,0) - y(0,0) - y(3,0);
1050 c(1,0) = x(1,0) - y(1,0) - y(4,0);
1051 c(2,0) = x(2,0) - y(2,0) - y(5,0);
1052 c(3,0) = pi0E - gammaE0 - gammaE1;
1053 if(_verbose>1) cerr << "set up the constraints " << endl;
1054
1055 if(_verbose>1)
1056 {
1057 cerr << "y: " << endl;
1058 for(int j=0;j<y.GetNrows();j++)
1059 {
1060 for(int k=0;k<y.GetNcols();k++)
1061 {
1062 cerr << j << " " << k << " " << y(j,k) << " ";
1063 }
1064 cerr << endl;
1065 }
1066 }
1067
1068 if(_verbose>1)
1069 {
1070 cerr << "x: " << endl;
1071 for(int j=0;j<x.GetNrows();j++)
1072 {
1073 for(int k=0;k<x.GetNcols();k++)
1074 {
1075 cerr << j << " " << k << " " << x(j,k) << " ";
1076 }
1077 cerr << endl;
1078 }
1079 }
1080
1081 if(_verbose>1)
1082 {
1083 cerr << "c: " << endl;
1084 for(int j=0;j<c.GetNrows();j++)
1085 {
1086 for(int k=0;k<c.GetNcols();k++)
1087 {
1088 cerr << j << " " << k << " " << c(j,k) << " ";
1089 }
1090 cerr << endl;
1091 }
1092 }
1093
1094
1095 //This is for 4 constraints
1096 // Derivatives of constraint equations wrt unknown quantities.
1097 // In this case the 4-momentum of the particle which decayed to the
1098 // two photons.
1099 a(0,0) = 1.0;
1100 a(1,1) = 1.0;
1101 a(2,2) = 1.0;
1102 for(i = 0; i < 3; i++)
1103 {
1104 a(3,i) = x(i,0)/pi0E;
1105 }
1106
1107 if(_verbose>1)
1108 {
1109 cerr << "a: " << endl;
1110 for(int j=0;j<a.GetNrows();j++)
1111 {
1112 for(int k=0;k<a.GetNcols();k++)
1113 {
1114 cerr << a(j,k) << " ";
1115 }
1116 cerr << endl;
1117 }
1118 }
1119
1120 aT.Transpose(a); // set a transpose
1121
1122 if(_verbose>1)
1123 {
1124 cerr << "aT: " << endl;
1125 for(int j=0;j<aT.GetNrows();j++)
1126 {
1127 for(int k=0;k<aT.GetNcols();k++)
1128 {
1129 cerr << aT(j,k) << " ";
1130 }
1131 cerr << endl;
1132 }
1133 }
1134
1135 // Set up derivatives matrix b(i,j) = dc(i)/dy(j):
1136 //This is for 4 constraints
1137 b(0,0) = -1.0;
1138 b(1,1) = -1.0;
1139 b(2,2) = -1.0;
1140 b(0,3) = -1.0;
1141 b(1,4) = -1.0;
1142 b(2,5) = -1.0;
1143 for(i = 0; i < 3; i++)
1144 {
1145 b(3,i) = -y(i,0)/gammaE0;
1146 b(3,i+3) = -y(i+3,0)/gammaE1;
1147 }
1148
1149 if(_verbose>1)
1150 {
1151 cerr << "b: " << endl;
1152 for(int j=0;j<b.GetNrows();j++)
1153 {
1154 for(int k=0;k<b.GetNcols();k++)
1155 {
1156 cerr << b(j,k) << " ";
1157 }
1158 cerr << endl;
1159 }
1160 }
1161
1162 bT.Transpose(b); // set b transpose
1163
1164 if(_verbose>1)
1165 {
1166 cerr << "bT: " << endl;
1167 for(int j=0;j<bT.GetNrows();j++)
1168 {
1169 for(int k=0;k<bT.GetNcols();k++)
1170 {
1171 cerr << bT(j,k) << " ";
1172 }
1173 cerr << endl;
1174 }
1175 }
1176
1177 // gb is a utility matrix we'll need a lot below
1178 gb = b * _cov * bT;
1179 if(_verbose>1)
1180 {
1181 cerr << "gb: " << endl;
1182 for(int j=0;j<gb.GetNrows();j++)
1183 {
1184 for(int k=0;k<gb.GetNcols();k++)
1185 {
1186 cerr << gb(j,k) << " ";
1187 }
1188 cerr << endl;
1189 }
1190 }
1191
1192 if( gb.Determinant() == 0 ) break; // for the case when gb is singular
1193
1194 gb.Invert();
1195 if(_verbose>1)
1196 {
1197 cerr << "gb inverted: " << endl;
1198 for(int j=0;j<gb.GetNrows();j++)
1199 {
1200 for(int k=0;k<gb.GetNcols();k++)
1201 {
1202 cerr << gb(j,k) << " ";
1203 }
1204 cerr << endl;
1205 }
1206
1207 cerr << "(aT * gb * a): " << endl;
1208 for(int j=0;j<(aT * gb * a).GetNrows();j++)
1209 {
1210 for(int k=0;k<(aT * gb * a).GetNcols();k++)
1211 {
1212 cerr << (aT * gb * a)(j,k) << " ";
1213 }
1214 cerr << endl;
1215 }
1216 }
1217
1218 // Is this what we should do?
1219 if( (aT * gb * a).Determinant() == 0.0 ) break;
1220 if(_verbose>1)
1221 {
1222 cerr << "Before invert: " << endl;
1223 cerr << "Determinant: " << (aT * gb * a).Determinant() << endl;
1224 (aT * gb * a).Invert();
1225 cerr << "After invert: " << endl;
1226 }
1227
1228 //if(_missingParticle)
1229 {
1230 // calculate the change to be made to the missing particle quantities
1231 xsi = ((aT * gb * a).Invert())*(aT * gb * c);
1232 // update x
1233 x -= xsi;
1234 // calculate the changes to be made to the measured quantities
1235 delta = _cov * bT * gb * (c - a*xsi);
1236
1237 if(_verbose>1)
1238 {
1239 cerr << "xsi: " << endl;
1240 for(int j=0;j<xsi.GetNrows();j++)
1241 {
1242 for(int k=0;k<xsi.GetNcols();k++)
1243 {
1244 cerr << xsi(j,k) << " ";
1245 }
1246 cerr << endl;
1247 }
1248
1249 cerr << "(c - a*xsi): " << endl;
1250 for(int j=0;j<(c - a*xsi).GetNrows();j++)
1251 {
1252 for(int k=0;k<(c - a*xsi).GetNcols();k++)
1253 {
1254 cerr << (c - a*xsi)(j,k) << " ";
1255 }
1256 cerr << endl;
1257 }
1258
1259 }
1260
1261 }
1262 //else delta = _cov * bT * gb * c;
1263
1264 // update the measured quantities y
1265 y -= delta;
1266
1267 if(_verbose>1)
1268 {
1269 cerr << "delta: " << endl;
1270 for(int j=0;j<delta.GetNrows();j++)
1271 {
1272 for(int k=0;k<delta.GetNcols();k++)
1273 {
1274 cerr << delta(j,k) << " ";
1275 }
1276 cerr << endl;
1277 }
1278 }
1279
1280 if(_verbose>1)
1281 {
1282 cerr << "new y: " << endl;
1283 for(int j=0;j<y.GetNrows();j++)
1284 {
1285 for(int k=0;k<y.GetNcols();k++)
1286 {
1287 cerr << y(j,k) << " ";
1288 }
1289 cerr << endl;
1290 }
1291 }
1292
1293
1294 eps = yi - y;
1295 epsT.Transpose(eps);
1296 // if(TMath::Prob((epsT * bT * gb * b * eps)(0,0),_ndf) < 1.e-10) break;
1297 if(_verbose>1)
1298 {
1299 cerr << "eps: " << endl;
1300 for(int j=0;j<eps.GetNrows();j++)
1301 {
1302 for(int k=0;k<eps.GetNcols();k++)
1303 {
1304 cerr << eps(j,k) << " ";
1305 }
1306 cerr << endl;
1307 }
1308 }
1309
1310
1311 } /* end iteration process */
1312
1313 // get the total changes in the measured quantities
1314 eps = yi - y;
1315 epsT.Transpose(eps);
1316
1317 // get the fit covariance matrix
1318 //if(_missingParticle)
1319 if( (aT * gb * a).Determinant() != 0.0 )
1320 {
1321 cov_fit = _cov - _cov*bT*gb*b*_cov
1322 + _cov*bT*gb*a*((aT*gb*a).Invert())*aT*gb*b*_cov;
1323 }
1324 else
1325 //cov_fit = _cov - _cov*bT*gb*b*_cov;
1326 {
1327 if(_verbose>1) cerr << "Not going to fill the missing cov matrix output 'cos is singular...." << endl;
1328 }
1329
1330 // get chi2 and the pulls
1331 _pulls.resize(dim);
1332 for(i = 0; i < dim; i++)
1333 {
1334 _pulls[i] = -eps(i,0)/sqrt(_cov(i,i) - cov_fit(i,i));
1335 }
1336
1337 _chi2 = (epsT * bT * gb * b * eps)(0,0); // the (0,0)...only...element
1338
1339 if(_verbose>0) cerr << "Prob: " << TMath::Prob(_chi2, _ndf) << " " << _chi2 << " " << _ndf << endl;
1340
1341 // update output of measured quantities
1342 if(_verbose>1) cerr << "Filling the 4vecs output...." << endl;
1343 for(i = 0; i < (int)_kDataFinal_out.size(); i++)
1344 {
1345 _kDataFinal_out[i].setMass(0.0);
1346 _kDataFinal_out[i].setMomentum(DVector3(y(3*i+0,0), y(3*i+1,0), y(3*i+2,0)));
1347 }
1348 if(_verbose>1) cerr << "Filled the 4vecs output...." << endl;
1349
1350 // set the missing particle covariance matrix
1351 if(_verbose>1) cerr << "Filling the missing cov matrix output...." << endl;
1352 // Not sure about this part
1353 if( (aT * gb * a).Determinant() != 0.0 )
1354 {
1355 if(_missingParticle)
1356 this->_SetMissingParticleErrors((aT * gb * a).Invert(),x);
1357 if(_verbose>1) cerr << "Filled the missing cov matrix output...." << endl;
1358 }
1359 else
1360 {
1361 if(_verbose>1) cerr << "Not going to fill the missing cov matrix output 'cos singular...." << endl;
1362 }
1363}
1364//_____________________________________________________________________________
1365
1366//_____________________________________________________________________________
1367
1368void DKinFit::_ResetForNewFit(){
1369 /// Reset the DKinFit object for a new fit (private function).
1370 // 1st make sure vectors are right size for current p4in
1371 _kDataInitial_in.clear();
1372 _kDataFinal_in.clear();
1373 _kDataInitial_out.clear();
1374 _kDataFinal_out.clear();
1375 _pulls.clear();
1376
1377 // reset bools to default values
1378 _missingParticle = false;
1379 _extraC_miss.clear();
1380}
1381//_____________________________________________________________________________
1382
1383double DKinFit::NameToMass(const string &__name) {
1384 /// Function used to translate names to masses for missing particles.
1385 /**
1386 * This function is used to translate @a name to a mass. It is used by
1387 * Fit to change a missing particle name to a missing mass. All names are
1388 * lower case and have no white space,-'s or _'s.
1389 *
1390 * The function first makes all letters lowercase, then removes any white
1391 * space, -'s or _'s from the name. It checks the name against the list of
1392 * known names. If it doesn't find a match, it outputs a warning message and
1393 * returns 0.
1394 *
1395 */
1396 if(__name == string()) return -1.;
1397
1398 string name = __name;
1399 // make all letters lower case
1400 std::transform(name.begin(),name.end(),name.begin(),(int(*)(int))std::tolower);
1401 // remove any white space, -'s or _'s
1402 for(string::size_type pos = 0; pos < name.size(); pos++){
1403 if(name[pos] == ' ' || name[pos] == '-' || name[pos] == '_'){
1404 name.erase(pos,1);
1405 pos--;
1406 }
1407 }
1408
1409 // mesons
1410 if(name == "pi" || name == "pi+" || name == "pi-") return 0.13957;
1411 if(name == "pi0") return 0.13498;
1412 if(name == "k" || name == "k+" || name == "k-") return 0.493677;
1413 if(name == "k0" || name == "kshort" || name == "klong") return 0.497672;
1414 if(name == "eta") return 0.54730;
1415 if(name == "rho" || name == "rho+" || name == "rho-" || name == "rho0")
1416 return 0.7693;
1417 if(name == "omega") return 0.78257;
1418 if(name == "eta'" || name == "etaprime") return 0.95778;
1419 if(name == "phi") return 1.019417;
1420
1421 // baryons
1422 if(name == "p" || name == "proton" || name == "antiproton" || name == "pbar"
1423 || name == "p-" || name == "p+") return 0.93827;
1424 if(name == "neutron" || name == "n") return 0.93957;
1425
1426 // if we get here it's not on the list
1427 std::cout << "Warning! <DKinFit::NameToMass> Unknown particle name: "
1428 << __name << ". Mass will be returned as 0." << std::endl;
1429 return 0.;
1430}
1431//_____________________________________________________________________________
1432
1433void DKinFit::_SetToBadFit()
1434{
1435 /// Set fit output to @a bad values.
1436 /** Sets pulls and \f$ \chi^2 \f$ to 666, sets degrees of freedom and output
1437 * kinematic quantites to 0.
1438 */
1439 for(int i = 0; i < (int)_pulls.size(); i++) _pulls[i] = 666.;
1440 _chi2 = -1.;
1441 _ndf = 0;
1442 for(int i = 0; i < (int)_kDataInitial_out.size(); i++)
1443 {
1444 _kDataInitial_out[i].setMomentum(DVector3(-666.,-666.,-666.));
1445 _kDataInitial_out[i].setMass(-666.);
1446 }
1447 for(int i = 0; i < (int)_kDataFinal_out.size(); i++)
1448 {
1449 _kDataFinal_out[i].setMomentum(DVector3(-666.,-666.,-666.));
1450 _kDataFinal_out[i].setMass(-666.);
1451 }
1452}
1453
1454//_____________________________________________________________________________
1455void DKinFit::_SetMissingParticleErrors(const DMatrix &__missingCov, const DMatrix &__x){
1456 /// Calculate the missing particle errors from the last fit
1457
1458 double p = sqrt(__x(0,0)*__x(0,0) + __x(1,0)*__x(1,0) + __x(2,0)*__x(2,0));
1459 DLorentzVector p4(__x(0,0),__x(1,0),__x(2,0),sqrt(p*p + pow(_missingMass,2)));
1460
1461 // kinematic quantities in tracking coordinates
1462 //double phi = PhiTrack(p4);
1463 //double lam = LambdaTrack(p4);
1464 //double alpha = AlphaTrack(p4);
1465 double phi = 1.0;
1466 double lam = 1.0;
1467 double alpha = 1.0;
1468
1469 // missing particle errors in lab coordinates
1470 double sigma_px_2 = __missingCov(0,0);
1471 double sigma_py_2 = __missingCov(1,1);
1472 double sigma_pz_2 = __missingCov(2,2);
1473
1474 // jacobian elements we need
1475 double dp_dpx = __x(0,0)/p;
1476 double dp_dpy = __x(1,0)/p;
1477 double dp_dpz = __x(2,0)/p;
1478
1479 double dlam_dpx = (-sin(alpha) - sin(lam)*dp_dpx)/(p*cos(lam));
1480 double dlam_dpy = (cos(alpha) - sin(lam)*dp_dpy)/(p*cos(lam));
1481 double dlam_dpz = -tan(lam)*dp_dpz/p;
1482
1483 double dphi_dpx = (sin(lam)*p*cos(phi)*dlam_dpx - dp_dpx*cos(lam)*cos(phi))
1484 /(-sin(phi)*p*cos(lam));
1485 double dphi_dpy = (sin(lam)*p*cos(phi)*dlam_dpy - dp_dpy*cos(lam)*cos(phi))
1486 /(-sin(phi)*p*cos(lam));
1487 double dphi_dpz = (1 + sin(lam)*p*cos(phi)*dlam_dpz
1488 - dp_dpz*cos(lam)*cos(phi))/(-sin(phi)*p*cos(lam));
1489
1490 // get error on p
1491 _sigma_missing[0] = sqrt(sigma_px_2*(dp_dpx*dp_dpx)
1492 + sigma_py_2*(dp_dpy*dp_dpy)
1493 + sigma_pz_2*(dp_dpz*dp_dpz));
1494
1495 // get error on lambda
1496 _sigma_missing[1] = sqrt(sigma_px_2*(dlam_dpx*dlam_dpx)
1497 + sigma_py_2*(dlam_dpy*dlam_dpy)
1498 + sigma_pz_2*(dlam_dpz*dlam_dpz));
1499
1500 // get error on phi
1501 _sigma_missing[2] = sqrt(sigma_px_2*(dphi_dpx*dphi_dpx)
1502 + sigma_py_2*(dphi_dpy*dphi_dpy)
1503 + sigma_pz_2*(dphi_dpz*dphi_dpz));
1504
1505}
1506//_____________________________________________________________________________
1507/// Set whether or not there's a missing particle
1508void DKinFit::SetMissingParticle(const double __missingMass)
1509{
1510 _missingParticle = true;
1511 _missingMass = __missingMass;
1512}
1513
1514//_____________________________________________________________________________
1515/// Set each extra mass constraint
1516void DKinFit::SetMassConstraint(const double __constraintMass, const vector<int> __constrainedParticles)
1517{
1518 _extraC++;
1519 _constraintMasses.push_back(__constraintMass);
1520 _constraintParticles.push_back(__constrainedParticles); /// Only final state particles for now
1521 if(_verbose)
1522 {
1523 for(int i=0;i<(int)_constraintMasses.size();i++)
1524 {
1525 cerr << "Pushing back constraint mass: " << _constraintMasses[i] << endl;
1526 }
1527 }
1528}
1529
1530//_____________________________________________________________________________
1531//
1532DLorentzVector DKinFit::MissingMomentum(vector<const DKinematicData*> initial, vector<const DKinematicData*> final)
1533{
1534 /// Just a function to return missing momentum
1535 if(_verbose>1) cerr << "Here in missing momentum...." << endl;
1536 DLorentzVector ret;
1537 float E = 0;
1538 float px = 0;
1539 float py = 0;
1540 float pz = 0;
1541 for(int i = 0; i < (int)initial.size(); i++)
1542 {
1543 E += initial[i]->energy();
1544 px += initial[i]->px();
1545 py += initial[i]->py();
1546 pz += initial[i]->pz();
1547 }
1548 for(int i = 0; i < (int)final.size(); i++)
1549 {
1550 E -= final[i]->energy();
1551 px -= final[i]->px();
1552 py -= final[i]->py();
1553 pz -= final[i]->pz();
1554 }
1555
1556 ret.SetXYZT(px, py, pz, E);
1557
1558 return ret;
1559}
1560
1561// Transform the 5x5 tracking error matrix into a 7x7 error matrix in cartesian
1562// coordinates
1563 DMatrixDSym DKinFit::Get7x7ErrorMatrix(const double mass,const double vec[5],
1564 const DMatrixDSym &C){
1565 DMatrixDSym C7x7(7);
1566 DMatrix J(7,5);
1567 DMatrixDSym C5x5(C);
1568
1569 // State vector
1570 double q_over_pt=vec[0];
1571 double phi=vec[1];
1572 double tanl=vec[2];
1573 double D=vec[3];
1574
1575 double cosl=cos(atan(tanl));
1576 double pt=1./fabs(q_over_pt);
1577 double p=pt/cosl;
1578 double p_sq=p*p;
1579 double E=sqrt(mass*mass+p_sq);
1580 double pt_sq=pt*pt;
1581 double cosphi=cos(phi);
1582 double sinphi=sin(phi);
1583 double q=(q_over_pt>0)?1.:-1.;
1584
1585 J(state_Px,state_q_over_pt)=-q*pt_sq*cosphi;
1586 J(state_Px,state_phi)=-pt*sinphi;
1587
1588 J(state_Py,state_q_over_pt)=-q*pt_sq*sinphi;
1589 J(state_Py,state_phi)=pt*cosphi;
1590
1591 J(state_Pz,state_q_over_pt)=-q*pt_sq*tanl;
1592 J(state_Pz,state_tanl)=pt;
1593
1594 J(state_E,state_q_over_pt)=-q*pt*p_sq/E;
1595 J(state_E,state_tanl)=pt_sq*tanl/E;
1596
1597 J(state_X,state_phi)=-D*cosphi;
1598 J(state_X,state_D)=-sinphi;
1599
1600 J(state_Y,state_phi)=-D*sinphi;
1601 J(state_Y,state_D)=cosphi;
1602
1603 J(state_Z,state_z)=1.;
1604
1605 // C'= JCJ^T
1606 C7x7=C5x5.Similarity(J);
1607
1608 return C7x7;
1609}