1 | #include <stdio.h> |
2 | #include <stdlib.h> |
3 | |
4 | #include <cctype> |
5 | #include <algorithm> |
6 | #include "PID/DKinFit.h" |
7 | using namespace std; |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | DKinFit::DKinFit(){ |
31 | |
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 | |
44 | void DKinFit::_Copy(const DKinFit &__kfit){ |
45 | |
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 | |
67 | void DKinFit::Fit() |
68 | { |
69 | |
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; |
77 | const int dim = numypar*numParts; |
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 | |
86 | _cov.ResizeTo(dim, dim); |
87 | |
88 | |
89 | |
90 | if(_cov.GetNrows() != dim || _cov.GetNcols() != dim) |
91 | { |
92 | |
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]); |
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 | |
115 | |
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 | |
132 | |
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]); |
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 | |
159 | double chi2=1e16,chi2old=0.; |
160 | |
161 | |
162 | |
163 | int numConstraints = 4; |
164 | if(_extraC) numConstraints += _extraC; |
165 | |
166 | |
167 | |
168 | |
169 | _ndf = 4; |
170 | if(_missingParticle) _ndf -= 3; |
171 | _ndf += _extraC; |
172 | |
173 | |
174 | |
175 | |
176 | if(_verbose>0) cerr << "numConstraints/ndf: " << numConstraints << " " << _ndf << endl; |
177 | |
178 | int i; |
179 | |
180 | vector<double> mass(numParts); |
181 | vector<double> energy(numParts); |
182 | double missingEnergy = 999.0; |
183 | vector<int> initialOrFinal(numParts); |
184 | |
185 | |
186 | |
187 | |
188 | DVector3 p3_miss; |
189 | |
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 | |
193 | |
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 | |
202 | |
203 | |
204 | |
205 | DLorentzVector missingMomentum; |
206 | |
207 | |
208 | |
209 | DMatrix yi(dim,1),y(dim,1); |
210 | |
211 | DMatrix x(3,1); |
212 | |
213 | |
214 | |
215 | DMatrix a(numConstraints,3),aT(3,numConstraints); |
216 | |
217 | |
218 | |
219 | DMatrix b(numConstraints,dim),bT(dim,numConstraints); |
220 | |
221 | |
222 | DMatrix c(numConstraints,1); |
223 | DMatrix c0(numConstraints,1); |
224 | |
225 | DMatrix eps(dim,1),epsT(1,dim); |
226 | |
227 | |
228 | DMatrix delta(dim,1); |
229 | |
230 | DMatrix xsi(3,1); |
231 | |
232 | |
233 | DMatrix gb(numConstraints,numConstraints); |
234 | DMatrix gb0(numConstraints,numConstraints); |
235 | |
236 | DMatrix cov_fit(dim,dim); |
237 | |
238 | DMatrix L(numConstraints,1); |
239 | DMatrix LT(1,numConstraints); |
240 | |
241 | |
242 | |
243 | |
244 | for(i = 0; i < numInitial; i++) |
245 | { |
246 | mass[i] = _kDataInitial_in[i].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(); |
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; |
261 | |
262 | |
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 | |
283 | for(int iter = 0; iter < 10; iter++) |
284 | { |
285 | |
286 | |
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 | |
322 | |
323 | |
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 | |
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); |
364 | |
365 | for(int k=0;k<numConstraintParticles;k++) |
366 | { |
367 | int partIndex = _constraintParticles[j][k]; |
368 | int Eindex = numInitial + partIndex; |
369 | |
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; |
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 | |
407 | |
408 | |
409 | |
410 | |
411 | |
412 | |
413 | |
414 | |
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 | |
431 | |
432 | if(_missingParticle) |
433 | { |
434 | for(int j=0;j<3;j++) |
435 | { |
436 | a(j,j) = -1.0; |
437 | a(3,j) = -x(j,0)/missingEnergy; |
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); |
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 | |
480 | |
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 | |
503 | for(int k=0;k<numConstraintParticles;k++) |
504 | { |
505 | int partIndex = _constraintParticles[j][k]; |
506 | int Eindex = numInitial + partIndex; |
507 | |
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); |
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 | |
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; |
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 | |
599 | if(_verbose>1) cerr << "After invert: " << endl; |
600 | |
601 | if( (aT * gb * a).Determinant() == 0.0 ) break; |
602 | |
603 | |
604 | xsi = ((aT * gb * a).Invert())*(aT * gb * c); |
605 | |
606 | x -= xsi; |
607 | |
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 | |
641 | |
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 | |
675 | L=gb*c; |
676 | LT.Transpose(L); |
677 | |
678 | |
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 | } |
698 | |
699 | |
700 | eps = yi - y; |
701 | epsT.Transpose(eps); |
702 | |
703 | |
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 | |
721 | |
722 | |
723 | |
724 | |
725 | |
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 | |
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 | |
760 | _chi2=chi2; |
761 | |
762 | |
763 | |
764 | |
765 | |
766 | |
767 | if(_verbose>0) cerr << "Prob: " << TMath::Prob(_chi2, _ndf) << " " << _chi2 << " " << _ndf << endl; |
768 | |
769 | |
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 | |
788 | |
789 | for(i = 0; i < numFinal; i++) |
790 | { |
791 | if(_verbose>=1) cerr << "detector: " << i << " " << _kDataFinal_in[i].t1_detector() << endl; |
792 | |
793 | |
794 | |
795 | if(_kDataFinal_in[i].t1_detector() == SYS_BCAL || _kDataFinal_in[i].t1_detector() == SYS_TOF) |
796 | { |
797 | |
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 | |
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 | |
824 | } |
825 | } |
826 | } |
827 | |
828 | if(_verbose>1) cerr << "Filling the missing cov matrix output...." << endl; |
829 | |
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 | |
846 | void DKinFit::FitTwoGammas(const float __missingMass, const float errmatrixweight=1.0) |
847 | { |
848 | |
849 | |
850 | |
851 | |
852 | |
853 | |
854 | |
855 | _missingMass = __missingMass; |
856 | _kDataInitial_out.clear(); |
857 | _kDataFinal_out.clear(); |
858 | |
859 | const int numInitial = (int)_kDataInitial_in.size(); |
860 | const int numFinal = (int)_kDataFinal_in.size(); |
861 | const int numParts = numInitial + numFinal; |
862 | |
863 | const int numypar = 3; |
864 | const int dim = numypar*numParts; |
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 | |
873 | _cov.ResizeTo(dim, dim); |
874 | |
875 | |
876 | |
877 | if(_cov.GetNrows() != dim || _cov.GetNcols() != dim) |
878 | { |
879 | |
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 | |
889 | { |
890 | for(int i=0;i<(int)_kDataInitial_in.size();i++) |
891 | { |
892 | _kDataInitial_out.push_back(_kDataInitial_in[i]); |
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]); |
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 | |
929 | |
930 | |
931 | |
932 | int numConstraints = 4; |
933 | if(_extraC) numConstraints = 5; |
934 | |
935 | numConstraints = 4; |
936 | |
937 | |
938 | _ndf = 4; |
939 | if(_missingParticle) _ndf -= 3; |
940 | if(_extraC) _ndf += 1; |
941 | |
942 | _ndf = 1; |
943 | |
944 | if(_verbose>0) cerr << "numConstraints/ndf: " << numConstraints << " " << _ndf << endl; |
945 | |
946 | int i; |
947 | |
948 | vector<double> mass(numParts); |
949 | vector<int> initialOrFinal(numParts); |
950 | |
951 | |
952 | |
953 | |
954 | |
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 | |
964 | |
965 | |
966 | |
967 | DLorentzVector missingMomentum; |
968 | |
969 | |
970 | |
971 | DMatrix yi(dim,1),y(dim,1); |
972 | |
973 | DMatrix x(3,1); |
974 | |
975 | |
976 | |
977 | DMatrix a(numConstraints,3),aT(3,numConstraints); |
978 | |
979 | |
980 | |
981 | DMatrix b(numConstraints,dim),bT(dim,numConstraints); |
982 | |
983 | |
984 | DMatrix c(numConstraints,1); |
985 | |
986 | DMatrix eps(dim,1),epsT(1,dim); |
987 | |
988 | |
989 | DMatrix delta(dim,1); |
990 | |
991 | DMatrix xsi(3,1); |
992 | |
993 | |
994 | DMatrix gb(numConstraints,numConstraints); |
995 | |
996 | DMatrix cov_fit(dim,dim); |
997 | |
998 | |
999 | |
1000 | |
1001 | for(i = 0; i < numInitial; i++) |
1002 | { |
1003 | mass[i] = _kDataInitial_in[i].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(); |
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; |
1016 | |
1017 | |
1018 | |
1019 | { |
1020 | |
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 | |
1027 | for(int iter = 0; iter < 10; iter++) |
1028 | { |
1029 | |
1030 | |
1031 | a.Zero(); |
1032 | aT.Zero(); |
1033 | b.Zero(); |
1034 | bT.Zero(); |
1035 | c.Zero(); |
1036 | delta.Zero(); |
1037 | xsi.Zero(); |
1038 | |
1039 | |
1040 | |
1041 | |
1042 | for(int j=0;j<numConstraints;j++) c(j,0) = 0.; |
1043 | |
1044 | |
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 | |
1096 | |
1097 | |
1098 | |
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); |
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 | |
1136 | |
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); |
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 | |
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; |
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 | |
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 | |
1229 | { |
1230 | |
1231 | xsi = ((aT * gb * a).Invert())*(aT * gb * c); |
1232 | |
1233 | x -= xsi; |
1234 | |
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 | |
1263 | |
1264 | |
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 | |
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 | } |
1312 | |
1313 | |
1314 | eps = yi - y; |
1315 | epsT.Transpose(eps); |
1316 | |
1317 | |
1318 | |
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 | |
1326 | { |
1327 | if(_verbose>1) cerr << "Not going to fill the missing cov matrix output 'cos is singular...." << endl; |
1328 | } |
1329 | |
1330 | |
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); |
1338 | |
1339 | if(_verbose>0) cerr << "Prob: " << TMath::Prob(_chi2, _ndf) << " " << _chi2 << " " << _ndf << endl; |
1340 | |
1341 | |
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 | |
1351 | if(_verbose>1) cerr << "Filling the missing cov matrix output...." << endl; |
1352 | |
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 | |
1368 | void DKinFit::_ResetForNewFit(){ |
1369 | |
1370 | |
1371 | _kDataInitial_in.clear(); |
1372 | _kDataFinal_in.clear(); |
1373 | _kDataInitial_out.clear(); |
1374 | _kDataFinal_out.clear(); |
1375 | _pulls.clear(); |
1376 | |
1377 | |
1378 | _missingParticle = false; |
1379 | _extraC_miss.clear(); |
1380 | } |
1381 | |
1382 | |
1383 | double DKinFit::NameToMass(const string &__name) { |
1384 | |
1385 | |
1386 | |
1387 | |
1388 | |
1389 | |
1390 | |
1391 | |
1392 | |
1393 | |
1394 | |
1395 | |
1396 | if(__name == string()) return -1.; |
1397 | |
1398 | string name = __name; |
1399 | |
1400 | std::transform(name.begin(),name.end(),name.begin(),(int(*)(int))std::tolower); |
1401 | |
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 | |
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 | |
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 | |
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 | |
1433 | void DKinFit::_SetToBadFit() |
1434 | { |
1435 | |
1436 | |
1437 | |
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 | |
1455 | void DKinFit::_SetMissingParticleErrors(const DMatrix &__missingCov, const DMatrix &__x){ |
1456 | |
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 | |
1462 | |
1463 | |
1464 | |
1465 | double phi = 1.0; |
1466 | double lam = 1.0; |
1467 | double alpha = 1.0; |
1468 | |
1469 | |
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 | |
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 | |
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 | |
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 | |
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 | |
1508 | void DKinFit::SetMissingParticle(const double __missingMass) |
1509 | { |
1510 | _missingParticle = true; |
1511 | _missingMass = __missingMass; |
1512 | } |
1513 | |
1514 | |
1515 | |
1516 | void DKinFit::SetMassConstraint(const double __constraintMass, const vector<int> __constrainedParticles) |
1517 | { |
1518 | _extraC++; |
1519 | _constraintMasses.push_back(__constraintMass); |
1520 | _constraintParticles.push_back(__constrainedParticles); |
1521 | if(_verbose) |
| |
1522 | { |
1523 | for(int i=0;i<(int)_constraintMasses.size();i++) |
| 2 | | Loop condition is true. Entering loop body | |
|
1524 | { |
1525 | cerr << "Pushing back constraint mass: " << _constraintMasses[i] << endl; |
| 3 | | Dereference of null pointer |
|
1526 | } |
1527 | } |
1528 | } |
1529 | |
1530 | |
1531 | |
1532 | DLorentzVector DKinFit::MissingMomentum(vector<const DKinematicData*> initial, vector<const DKinematicData*> final) |
1533 | { |
1534 | |
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 | |
1562 | |
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 | |
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 | |
1606 | C7x7=C5x5.Similarity(J); |
1607 | |
1608 | return C7x7; |
1609 | } |