215 {
216
217
218
219
222
223
225
226
228 unsigned nCores = cores.length();
230
231
232
233 bool flag2D = false;
234 if ((nStereoCores == 0) && (nCores > 3)) flag2D = true;
235 else if ((nStereoCores < 2) || (nCores - nStereoCores < 3))
237
238
241
242
252 double chi2;
254 double factor = 1.0;
255 int err = 0;
258
259
260 unsigned nTrial = 0;
261 while(nTrial < 100){
262
263
264 chi2 = 0;
265 for (unsigned j=0;j<4;j++) dchi2da[j]=0;
266 d2chi2d2a=zero4;
267
268
269 unsigned i=0;
270 while(
TMLink* l = cores[i++]){
272
273
275 const HepPoint3D& onTrack=l->positionOnTrack();
278 if (onWire.cross(onTrack).z() < 0.) leftRight =
WireHitLeft;
279
280
281 double distance;
282 double eDistance;
283 drift(
t, * l, t0Offset, distance, eDistance);
284 l->drift(distance,0);
285 l->drift(distance,1);
286 l->dDrift(eDistance,0);
287 l->dDrift(eDistance,1);
288 double eDistance2 = eDistance * eDistance;
289
290
291
293 double vmag =
v.mag();
294 double dDistance = vmag - distance;
295
297
298 this->dxda(*l,
t, dxda, dyda, dzda, vw);
299
300
301 dDda = (vmag > 0.)
302 ? ((
v.x() * (1. - vw.x() * vw.x()) -
303 v.y() * vw.x() * vw.y() -
v.z() * vw.x() * vw.z())
304 * dxda +
305 (
v.y() * (1. - vw.y() * vw.y()) -
306 v.z() * vw.y() * vw.z() -
v.x() * vw.y() * vw.x())
307 * dyda +
308 (
v.z() * (1. - vw.z() * vw.z()) -
309 v.x() * vw.z() * vw.x() -
v.y() * vw.z() * vw.y())
310 * dzda) / vmag :
Vector(4,0);
311 if (vmag <= 0.0) {
312 std::cout << " in fit " << onTrack << ", " << onWire;
314 }
315 dchi2da += (dDistance / eDistance2) * dDda;
316 d2chi2d2a += vT_times_v(dDda) / eDistance2;
317 double pChi2 = dDistance * dDistance / eDistance2;
318 chi2 += pChi2;
319
320
321 l->update(onTrack, onWire, leftRight, pChi2);
322 }
323
324
325 double change = chi2Old - chi2;
326
327 if (fabs(change) < 1.0e-5) break;
328 if (change < 0.) {
329 a += factor * da;
330 factor *= 0.1;
331 }else{
332 chi2Old = chi2;
333 if(flag2D){
334 f = dchi2da.sub(1,2);
335 e = d2chi2d2a.sub(1,2);
336 f = solve(e,f);
337 da[0]=f[0];
338 da[1]=f[1];
339 da[2]= 0;
340 da[3]= 0;
341 }else{
342
343 da = solve(d2chi2d2a, dchi2da);
344 }
345 }
346 a -= factor * da;
348 ++nTrial;
349 }
350
351
353 unsigned dim;
354 if(flag2D){
355 dim=2;
357 Eb = d2chi2d2a.sub(1,3);
358 Ec = Eb.inverse(err);
359 Ea[0][0] = Ec[0][0];
360 Ea[0][1] = Ec[0][1];
361 Ea[0][2] = Ec[0][2];
362 Ea[1][1] = Ec[1][1];
363 Ea[1][2] = Ec[1][2];
364 Ea[2][2] = Ec[2][2];
365 }else{
366 dim=4;
367 Ea = d2chi2d2a.inverse(err);
368 }
369
370
371 if(! err){
376 }else{
378 }
379
380 t._ndf = nCores - dim;
382
383
385
386 return err;
387}
CLHEP::HepSymMatrix SymMatrix
const HepPoint3D ORIGIN
Constants.
#define TFitAlreadyFitted
unsigned NStereoHits(const AList< TMLink > &links)
returns # of stereo hits.
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
A class to represent a track in tracking.
void dump(const std::string &message=std::string(""), const std::string &prefix=std::string("")) const
dumps debug information.
A class to relate TMDCWireHit and TTrack objects.
virtual unsigned objectType(void) const
returns object type.