StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TCFit.cxx
1 // @(#)root/base:$Name: $:$Id: TCFit.cxx,v 1.6 2010/01/28 18:19:05 perev Exp $
2 // Author: Victor Perev 05/08/03
3 
4 
5 #include <stdio.h>
6 #include <stdlib.h>
7 #include <assert.h>
8 #include <math.h>
9 #include "TError.h"
10 #include "TCFit.h"
11 #include "TNumDeriv.h"
12 #include "TVector3.h"
13 #include "TRandom.h"
14 #include "TLorentzVector.h"
15 #include "THelixTrack.h"
16 #if ROOT_VERSION_CODE < 331013
17 #include "TCL.h"
18 #else
19 #include "TCernLib.h"
20 #endif
21 ClassImp(TCFit)
22 //______________________________________________________________________________
23 TCFit::TCFit(const char *name,TCFitData *dat) :TNamed(name,"")
24 {
25  fDebug = 1;
26  SetData(dat);
27  for (int i=0;i<5;i++) (&fBigM)[i]=new TMatrixD;
28 
29 }
30 //______________________________________________________________________________
31 void TCFit::Reset()
32 {
33  for (int i=0;i<5;i++) {*((&fBigM)[i])=0.;}
34  fIter =0;
35  fCuts =0;
36  fMaxIter =500;
37  fMaxCuts =6;
38  fUPars = 0;
39  fUMeas = 0;
40  fUCons = 0;
41 }
42 //______________________________________________________________________________
43 TCFit::~TCFit()
44 {
45  for (int i=0;i<4;i++) delete (&fBigM)[i];
46 }
47 //______________________________________________________________________________
48 int TCFit::SetData(TCFitData *dat)
49 {
50  fDat = dat;
51  fDat->SetFitter(this);
52  return 0;
53 }
54 //______________________________________________________________________________
55 int TCFit::CheckOut()
56 {
57  int nBig = fUPars+fUCons;
58 
59  fDat->Evaluate();
60 
61  if (!fIter) {
62  fBigM->ResizeTo(nBig,nBig); (*fBigM)*=0.0;
63  fBigI->ResizeTo(nBig,nBig);
64  fBigB->ResizeTo(nBig,1);
65  fOldP->ResizeTo(nBig,1);
66  fAddP->ResizeTo(nBig,1);
67  }
68  for (int j1=0;j1<fUMeas; j1++) {
69  int i1 = fDat->GetId(j1);
70  if (!fIter) (*fOldP )[j1][ 0] = fDat->GetPar(i1);
71  (*fBigB )[j1][ 0] = -fDat->DFcn(i1);
72  (*fBigM)[j1][j1] = fDat->DDFcn(i1,i1);
73  for (int j2=0;j2<j1; j2++) {
74  int i2 = fDat->GetId(j2);
75  (*fBigM)[j1][j2] = fDat->DDFcn(i1,i2);
76  (*fBigM)[j2][j1] = (*fBigM)[j1][j2];
77  }
78  }
79  for (int jc=0;jc<fUCons; jc++) {
80  int ic = fDat->GetId(jc+fUPars);
81  double cc = fDat->GetPar(ic);
82  (*fBigB )[jc+fUPars][0] = -cc;
83  for (int jx=0;jx<fUPars; jx++) {
84  int ix = fDat->GetId(jx);
85  double cx = fDat->DCon(ic,ix);
86  (*fBigM)[fUPars+jc][jx] = cx;
87  (*fBigM)[jx][fUPars+jc] = cx;
88  } }
89 
90 
91  for (int jx=0; jx<fUPars;jx++) {
92  for (int jc=0;jc<fUCons;jc++) {
93  (*fBigB )[jx][ 0] -= (*fOldP)[fUPars+jc][0]*(*fBigM)[jx][fUPars+jc];
94  } }
95 
96  if (fDebug>=3) {
97  fBigM->Print("CheckOut.BigM");
98  fBigB->Print("CheckOut.BigB");
99  fDat->Print("CheckOut.dat");
100  }
101 
102 
103  return 0;
104 }
105 //______________________________________________________________________________
106 int TCFit::CheckIn()
107 {
108  for (int jx=0;jx<fUPars;jx++) {
109  int ix = fDat->GetId(jx);
110  fDat->GetPar(ix) = (*fOldP)[jx][0] + (*fAddP)[jx][0];
111  }
112  fDat->Evaluate();
113 // fDat->Print("CheckIn");
114  return 0;
115 }
116 //______________________________________________________________________________
117 int TCFit::Fit()
118 {
119  if (fDat->Ready()) return kBadAprx;
120  fUPars = fDat->GetUPars();
121  fUMeas = fDat->GetUMeas();
122  fUCons = fDat->GetUCons();
123 
124  fCuts = 0; fFitRes=0;
125 
126  fFcnQA[0]=3e33; fFcnQA[1]=3e33;
127  fConQA[0]=3e33; fConQA[1]=3e33;
128  fAddQA[0]=3e33; fAddQA[1]=3e33;
129 
130 
131 
132 
133  for (fIter=0; fIter < fMaxIter; fIter++)
134  {
135  fAkt = CheckStep();
136  switch (fAkt) {
137 
138  case kNextStep: (*fOldP)+=(*fAddP); FitStep(); break;
139 
140  case kNextCut: CutStep(); break;
141 
142  case kEndFit: EndStep(); return fFitRes;
143 
144  case kBadFit: EndStep(); return fFitRes;
145 
146  default: assert(0);
147  }//end case
148  if (fDebug>=2) PriStep("OneIter");
149 
150  }
151  fFitRes=kTooIter;
152  EndStep();
153  return fFitRes;
154 }
155 //______________________________________________________________________________
156 int TCFit::FitStep()
157 {
158  CheckOut();
159  fFcnQA[1]=fFcnQA[0];
160  fConQA[1]=fConQA[0];
161  fAddQA[1]=fAddQA[0];
162 
163  (*fBigI) = (*fBigM);
164  double det=0;
165  (*fBigI).Invert(&det);
166  (*fAddP) = (*fBigI)*(*fBigB);
167 // fAddP->Print("AddP");
168  return 0;
169 }
170 //______________________________________________________________________________
171 int TCFit::CutStep()
172 {
173  fCuts++; fIter--;
174  (*fAddP)*=0.5; return 0;
175 }
176 //______________________________________________________________________________
177 int TCFit::EndStep()
178 {
179  fDat->SetFail(fFitRes);
180  if (fDebug>=1) PriStep("EndStep");
181  if (fDebug>=2) fDat->Print("EndStep");
182  return 0;
183 }
184 //______________________________________________________________________________
185 int TCFit::CheckStep()
186 {
187  fFitRes = 0;
188 
189  if (fIter==0) return kNextStep;
190  CheckIn();
191  fFcnQA[0] = fDat->GetFcn();
192  fConQA[0] = 0;
193  fAddQA[0] = 0;
194  for (int jc=0;jc<fUCons; jc++) {
195  int ic = fDat->GetId(jc+fUPars);
196  fConQA[0]+= fabs(fDat->GetPar(ic)/fDat->GetTiny(ic));
197  }
198  if (fUCons>1) fConQA[0]/=fUCons;
199  for (int jx=0;jx<fUMeas; jx++) {
200  int ix = fDat->GetId(jx);
201  fAddQA[0]+= fabs((*fAddP)[jx][0]/fDat->GetTiny(ix));
202  }
203  fAddQA[0]/=fUMeas;
204 
205  if (fAddQA[0]<1. && fConQA[0] <1.) return kEndFit;
206 
207  if (fFcnQA[0] > fDat->GetBigFcn() ) fFitRes |=kBadFcn;
208  if (fConQA[0] > 1e10 ) fFitRes |=kBadCon;
209  if (fFitRes) return kBadFit;
210 
211  if (fAddQA[0]<fAddQA[1]
212  || fConQA[0]<fConQA[1]
213  || fFcnQA[0]<fFcnQA[1]
214  || fCuts>fMaxCuts) {
215  fCuts = 0;
216  fAddQA[1]=fAddQA[0];
217  fConQA[1]=fConQA[0];
218  fFcnQA[1]=fFcnQA[0];
219  return kNextStep;
220  }
221  return kNextCut;
222 }
223 //______________________________________________________________________________
224 double TCFit::ErMx(int jcol,int jrow) const
225 {
226  return (*fBigI)[jcol][jrow];
227 }
228 //______________________________________________________________________________
229 int TCFit::PriStep(const char *where)
230 {
231  if (where==0)where="";
232  printf("PriStep(%s) Iter=%d Cut=%d Akt=%d Fcn=%g(%g) Con=%g(%g) Add=%g(%g)\n"
233  ,where,fIter,fCuts,fAkt
234  ,fFcnQA[0],fFcnQA[1],fConQA[0],fConQA[1],fAddQA[0],fAddQA[1]);
235 
236  return 0;
237 }
238 
239 
240 
241 
242 
243 //______________________________________________________________________________
244 //______________________________________________________________________________
245 class TEST0 : public TCFitData {
246 public:
247  TEST0();
248 
249  double Fcn();
250  double Con(int icon);
251 
252  double MyDFcn(int ix);
253  double MyDDFcn(int ix);
254  void Update() {}
255  int Approx() {return 0;}
256 public:
257  double mX[3],mC;
258 
259 };
260 //______________________________________________________________________________
261 TEST0::TEST0():TCFitData("TEST0")
262 {
263  mX[0]=mX[1]=mX[2]=1.1;
264  mC=0;
265  AddPar(0,0, mX,3,"X",0.0001);
266  AddPar(2,3,&mC,1,"C",0.0001);
267 }
268 //______________________________________________________________________________
269 double TEST0::Fcn()
270 {
271  return mX[0]*mX[0]/1. + mX[1]*mX[1]/4 + mX[2]*mX[2]/9;
272 }
273 //______________________________________________________________________________
274 double TEST0::MyDFcn(int ix)
275 {
276  switch (ix) {
277  case 0: return 2*mX[0];
278  case 1: return 2*mX[1]/4;
279  case 2: return 2*mX[2]/9;
280  }
281  assert(0); return 0.;
282 
283 }
284 //______________________________________________________________________________
285 double TEST0::MyDDFcn(int ix)
286 {
287  switch (ix) {
288  case 0: return 2.;
289  case 1: return 2./4;
290  case 2: return 2./9;
291  }
292  assert(0); return 0.;
293 
294 }
295 //______________________________________________________________________________
296 double TEST0::Con(int)
297 {
298  return mX[2]-1.;
299 }
300 
301 
302 //______________________________________________________________________________
303 void TCFit::Test0()
304 {
305  TEST0 t0;
306  TCFit f0("Test0",&t0);
307  for (int itst=0;itst<3;itst++) {
308  double d1 = t0.DFcn(itst);
309  double d2 = t0.MyDFcn(itst);
310  double dd1 = t0.DDFcn(itst,itst);
311  double dd2 = t0.MyDDFcn(itst);
312  printf("%d %g %g %g %g\n",itst,d1,d2,dd1,dd2);
313  }
314  f0.Fit();
315 }
316 
317 
318 
319 ClassImp(TCFitData)
320 
321 class Deriv1st : public TNumDeriv {
322 public:
323  Deriv1st(TCFitData* dat);
324  Double_t Fcn(Double_t arg);
325  void SetKind(int kind) { fKind=kind;}
326 private:
327 TCFitData *fFitData;
328 int fKind;
329 };
330 //______________________________________________________________________________
331 Deriv1st::Deriv1st(TCFitData* dat):TNumDeriv("Deriv1st")
332 {
333  fFitData=dat;
334  fKind = -1;
335 }
336 //______________________________________________________________________________
337 Double_t Deriv1st::Fcn(Double_t arg)
338 {
339  double save = fFitData->GetPar(fIArg);
340  SetStep(fFitData->GetTiny(fIArg));
341  fFitData->GetPar(fIArg)+=arg;
342  fFitData->Update(); fFitData->Modify(0);
343  double ans = (fKind==(-1)) ? fFitData->Fcn() :fFitData->Con(fKind);
344  fFitData->GetPar(fIArg)=save;
345  fFitData->Update(); fFitData->Modify(0);
346  return ans;
347 }
348 
349 //______________________________________________________________________________
350 //______________________________________________________________________________
351 class Deriv2nd : public TNumDeriv {
352 public:
353  Deriv2nd(TCFitData* dat);
354  ~Deriv2nd();
355  Double_t Fcn(Double_t arg);
356  void SetKind(int kind) { fD1->SetKind(kind);}
357  void SetIJArg(int i,int j);
358 private:
359 TCFitData *fFitData;
360 Deriv1st *fD1;
361 };
362 //______________________________________________________________________________
363 Deriv2nd::Deriv2nd(TCFitData* dat):TNumDeriv("Deriv2nd")
364 {
365  fFitData = dat;
366  fD1 = new Deriv1st(fFitData);
367 }
368 //______________________________________________________________________________
369 Deriv2nd::~Deriv2nd()
370 {
371  delete fD1;
372 }
373 //______________________________________________________________________________
374 void Deriv2nd::SetIJArg(int i,int j)
375 {
376  this->SetIArg(i);
377  fD1 ->SetIArg(j);
378 }
379 //______________________________________________________________________________
380 Double_t Deriv2nd::Fcn(Double_t arg)
381 {
382  double save = fFitData->GetPar(fIArg);
383  fFitData->GetPar(fIArg)+=arg;
384  fFitData->Update(); fFitData->Modify(0);
385  double ans = fD1->DFcn(0.);
386  fFitData->GetPar(fIArg)=save;
387  fFitData->Update(); fFitData->Modify(0);
388  return ans;
389 }
390 //______________________________________________________________________________
391 //______________________________________________________________________________
392 TCFitData::TCFitData(const char *name, const char *title):TNamed(name,title)
393 {
394 
395  Reset();
396  fD1st = new Deriv1st(this);
397  fD2nd = new Deriv2nd(this);
398 }
399 //______________________________________________________________________________
400 void TCFitData::Reset()
401 {
402  memset(fBeg,0,fEnd-fBeg+1);
403  memset(fFixs,1,sizeof(fFixs));
404  for (int i=0; i<kMaxId;i++) fNams[i]="";
405  fFcn[1]=1e-6;
406  fFcn[2]=1e+6;
407 }
408 
409 //______________________________________________________________________________
410 TCFitData::~TCFitData()
411 {
412  delete fD1st;
413  delete fD2nd;
414 }
415 //______________________________________________________________________________
416 int TCFitData::AddPar(int tyPar, int idPar,double *par,int nPars,const char *name,double tiny)
417 {
418  assert(0<=tyPar && tyPar<=2);
419  assert(0<=idPar && idPar+nPars<kMaxId);
420  if (fFail ) return fFail;
421  if (!name) name="";
422  if (tiny<=0) tiny = 1e-3;
423 
424  for (int j=0;j<nPars;j++) {
425  TString ts;
426  if (name[0]) ts=name; else ts=idPar;
427  if (nPars>1) { ts +="["; ts +=j; ts+="]";}
428  fNams[idPar+j]= ts;
429  fTiny[idPar+j]= tiny;
430  fPars[idPar+j]= par+j;
431  fFixs[idPar+j]= 0;
432  fTyps[idPar+j]= tyPar;
433  }
434  fNPars[tyPar]+=nPars;
435  assert(fNPars[tyPar]<=kMaxId);
436  return 0;
437 }
438 //______________________________________________________________________________
439 Int_t TCFitData::GetId(const char *name) const //get par index by name
440 {
441  int i=0;
442  int n = fNPars[0]+fNPars[1]+fNPars[2];
443  for (int j=0;i<n;j++) {
444  if (!fNams[j].Length()) continue;
445  i++;
446  if (fNams[j]==name) return j;
447  }
448  printf("TCFitData::GetId(\"%s\") UNKNOWN name\n",name);
449  return -1;
450 }
451 //______________________________________________________________________________
452 const char *TCFitData::GetNam(Int_t idx) const //get par name by index
453 {
454  return fNams[idx].Data();
455 }
456 //______________________________________________________________________________
457 int TCFitData::GetType(int id) const
458 {
459 //get type(Meas,Slac,Constr) by id
460  return fTyps[id];
461 }
462 
463 
464 //______________________________________________________________________________
465 int TCFitData::GetId(int jdx) const
466 {
467  return fIndx[jdx];
468 }
469 //______________________________________________________________________________
470 int TCFitData::GetJd(int idx) const
471 {
472  return fJndx[idx];
473 }
474 
475 //______________________________________________________________________________
476 void TCFitData::FixPar(int idx,int yes) //(dis/en)able parameter
477 {
478  yes = (yes!=0);
479  assert(fNams[idx].Length());
480  if (fFixs[idx]==yes) return;
481  fFixs[idx]=yes;
482  fNFixs[fTyps[idx]]+= ( (yes)? 1:-1 );
483 
484 }
485 //______________________________________________________________________________
486 int TCFitData::IsFixed(int idx) const
487 { return fFixs[idx];}
488 
489 //______________________________________________________________________________
490 double TCFitData::GetPar (int idx) const
491 {return *fPars[idx];}
492 //______________________________________________________________________________
493 double &TCFitData::GetPar (int idx)
494 {
495  fModi=1; return *fPars[idx];
496 }
497 //______________________________________________________________________________
498 void TCFitData::Evaluate()
499 {
500  if (Modified()) Update(); Modify(0);
501  fFcn[0] = Fcn();
502  for (int jc=GetUPars();jc<GetUPars()+GetUCons();jc++) {
503  int ic = GetId(jc);
504  GetPar(ic) = Con(ic);
505  }
506 }
507 //______________________________________________________________________________
508 int TCFitData::Ready()
509 {
510  int n = 0;
511  for (int ity=0;ity<3;ity++) {
512  for (int i=0;i<kMaxId;i++) {
513  if (!fNams[i].Length()) continue;
514  if ( fFixs[i]) continue;
515  if ( fTyps[i]!=ity) continue;
516  fIndx[n]=i;
517  fJndx[i]=n++;
518  } }
519  assert(n==fNPars[0]-fNFixs[0]+fNPars[1]-fNFixs[1]+fNPars[2]-fNFixs[2]);
520  return Approx();
521 }
522 
523 //______________________________________________________________________________
524 Double_t TCFitData::DFcn(int ipar)
525 {
526  fD1st->SetKind(-1);
527  fD1st->SetIArg(ipar);
528  return fD1st->DFcn();
529 }
530 
531 //______________________________________________________________________________
532 Double_t TCFitData::DDFcn(int ipar,int jpar)
533 {
534  fD2nd->SetKind(-1);
535  fD2nd->SetIJArg(ipar,jpar);
536  return fD2nd->DFcn();
537 }
538 
539 //______________________________________________________________________________
540 Double_t TCFitData::Con(int)
541 {
542 assert(0);
543 }
544 
545 //______________________________________________________________________________
546 Double_t TCFitData::DCon(int icon,int ipar)
547 {
548  fD1st->SetKind(icon);
549  fD1st->SetIArg(ipar);
550  return fD1st->DFcn();
551 }
552 //______________________________________________________________________________
553 double TCFitData::ErMx(int icol,int irow) const
554 {
555  int jcol = GetJd(icol);
556  int jrow = GetJd(irow);
557  return fFitter->ErMx(jcol,jrow);
558 }
559 //______________________________________________________________________________
560 void TCFitData::Print(const char *name) const
561 {
562 static const char *ty[3]={"Meas","Slack","Cnsr"};
563 static const char *fx[2]={"Used","Fixed"};
564 
565  if (!name) name = "";
566  printf("TCFitData::Print(%s) nMeas=%d(%d) nSlac=%d(%d) nCons=%d(%d)\n"
567  ,name
568  ,GetNMeas(),GetUMeas()
569  ,GetNSlac(),GetUSlac()
570  ,GetNCons(),GetUCons());
571  for (int i=0;i<kMaxId;i++) {
572  if (!fNams[i].Length()) continue;
573  printf("%2d - %s\t",i,fNams[i].Data());
574  printf(" %s.%s ",ty[fTyps[i]],fx[fFixs[i]]);
575  printf(" %g \n",*fPars[i]);
576  }
577 
578 }
579 
580 
581 //______________________________________________________________________________
582 void TkPars::Print(const char *tit) const
583 {
584  if (!tit) tit="";
585  printf("TkPars::Print(%s) \tDca=%g Z=%g Phi=%g ptin=%g tanl=%g curv=%g"
586  ,tit,dca,z,phi,ptin,tanl,curv);
587 
588 }
589 //______________________________________________________________________________
590 void TkPars::Reset()
591 {
592  double tmp = hz;
593  memset(this,0,sizeof(TkPars));
594  hz = tmp;
595 }
596 //______________________________________________________________________________
597 void TkPars::SetHz(double fact)
598 {
599  hz = (0.000299792458 * 4.98478)*fact;
600 }
601 //______________________________________________________________________________
602 TkPars &TkPars::operator+=(const TkPars &a)
603 {
604  for (int i=0;i<5;i++) Arr()[i]+=a.Arr()[i];
605  if (phi<-M_PI) phi+=M_PI*2;
606  if (phi> M_PI) phi-=M_PI*2;
607  return *this;
608 }
609 //______________________________________________________________________________
610 TLorentzVector TkPars::P4() const
611 {
612  double pt = fabs(1./ptin);
613  double e = sqrt(pt*pt*(1.+tanl*tanl)+mass*mass);
614  TLorentzVector v(pt*cos(phi),pt*sin(phi),pt*tanl,e);
615  return v;
616 }
617 //______________________________________________________________________________
618 void TkPars::P4D(double D[4][5]) const
619 {
620  enum {Kdca,Kz,Kphi,Kptin,Ktanl};
621 
622  memset(D[0],0,4*5*sizeof(D[0][0]));
623  double pt = fabs(1./ptin);
624  double dpt = -pt/ptin;
625  double e = sqrt(pt*pt*(1.+tanl*tanl)+mass*mass);
626 // TLorentzVector v(pt*cos(phi),pt*sin(phi),pt*tanl,e);
627 
628  D[3][Kptin] = pt*(1.+tanl*tanl)/e*dpt;
629  D[3][Ktanl] = pt*pt*tanl/e*dpt;
630 
631  D[0][Kptin] = cos(phi)*dpt;
632  D[1][Kptin] = sin(phi)*dpt;
633  D[2][Kptin] = tanl*dpt;
634  D[2][Ktanl] = pt;
635 
636 }
637 //______________________________________________________________________________
638 void TkPars::Set(const TVector3 &v3,const TVector3 &d3,double pti)
639 {
640  dca = v3.Perp();
641  phi = d3.Phi();
642  z = v3[2];
643  if ((v3.Cross(d3)).Z()<0) dca=-dca;
644  ptin = pti;
645  tanl = d3.Pz()/d3.Pt();
646 }
647 //______________________________________________________________________________
648 //______________________________________________________________________________
649 void TkPars::Get(TVector3 *v3,TVector3 *d3,double *pts) const
650 {
651  if (v3 ) { v3->SetXYZ(dca*sin(phi),-dca*cos(phi),z) ;}
652  if (d3 ) { v3->SetXYZ(cos(phi),sin(phi),tanl); v3->SetMag(1.);}
653  if (pts) {*pts = 1./ptin ;}
654 }
655 //______________________________________________________________________________
656 TVector3 TkPars::V3() const
657 {
658  TVector3 v3; Get(&v3); return v3;
659 }
660 //______________________________________________________________________________
661 void TkPars::Fill(THelixTrack &hlx)
662 {
663  TLorentzVector p4 = P4();
664  TVector3 v3 = V3();
665  double h[6];
666  v3.GetXYZ(h);
667  p4.Vect().GetXYZ(h+3);
668  hlx.Set(h,h+3,curv);
669 
670 }
671 //______________________________________________________________________________
672 const char* TkPars::Name(int mem)
673 {
674 static const char* nams[]={ "Dca", "Z", "Phi", "Ptin", "TanL"};
675  return nams[mem];
676 }
677 //______________________________________________________________________________
678 double TkPars::Tiny(int mem)
679 {
680 static const double tiny[]={3.14/180/10., 0.01, 3.14/180/10, 0.01, 0.001};
681  return tiny[mem];
682 }
683 //______________________________________________________________________________
684 void TkPars::Rand(const TkErrs &errs)
685 {
686  for (int iv=0;iv<5;iv++) {
687  Arr()[iv]+= gRandom->Gaus(0.,sqrt(errs.Get(iv,iv)));
688  }
689  Update();
690 }
691 
692 //______________________________________________________________________________
693 //______________________________________________________________________________
694 void TkErrs::Reset()
695 {
696 static const double DX=0.3,DZ=0.3,DPhi=0.03,DPti=1.,DTan=0.05;
697  memset(emx,0,sizeof(emx));
698  double fak = 0.3*0.3;
699  emx[ 0]=DX*DX *fak;
700  emx[ 2]=DZ*DZ *fak;
701  emx[ 5]=DPhi*DPhi *fak;
702  emx[ 9]=DPti*DPti *fak;
703  emx[14]=DTan*DTan *fak;
704 }
705 //______________________________________________________________________________
706 void TkErrs::Set(int i,int j,double err)
707 {
708  if (i<j) { int ii=i; i=j; j=ii;}
709  emx[((i+1)*i)/2+j]=err;
710 }
711 //______________________________________________________________________________
712 double TkErrs::Get(int i,int j) const
713 {
714  if (i<j) { int ii=i; i=j; j=ii;}
715  return emx[((i+1)*i)/2+j];
716 }
717 //______________________________________________________________________________
718 void TkErrs::Invert()
719 {
720  TCL::trsinv(emx,emx,5);
721 }
722 //______________________________________________________________________________
723 double TkErrs::Xi2(const TkPars &pars) const
724 {
725  double chi2;
726  TCL::trasat(pars.Arr(),emx,&chi2,1,5);
727  return chi2;
728 }
729 //______________________________________________________________________________
730 void TkErrs::Mpy(const TkPars &pars,double der[5]) const
731 {
732  TCL::trsa(emx, pars.Arr() ,der,5,1);
733 }
734 //______________________________________________________________________________
735 //______________________________________________________________________________
736 TVector3 VxPars::V3() const
737 { return TVector3(x); }
738 
739 //______________________________________________________________________________
740 //______________________________________________________________________________
741 ClassImp(TCFitV0)
742 //______________________________________________________________________________
743 TCFitV0::TCFitV0():TCFitData("TCFitV0","Two prong decay fit")
744 {
745  Reset();
746 }
747 //______________________________________________________________________________
748 void TCFitV0::Reset()
749 {
750  for (int i=0;i<2;i++) {
751  mTkBas[i].Reset();
752  mTEBas[i].Reset();
753  mTkFit[i].Reset();
754  mTkDif[i].Reset();
755  }
756  memset(mBeg,0,mEnd-mBeg+1);
757  TCFitData::Reset();
758 }
759 //______________________________________________________________________________
760 int TCFitV0::Ready()
761 {
762  if (!mReady) {
763  mReady=46;
764  TCFitData::Reset();
765 
766  for (int itk=0;itk<2;itk++) {
767  mTkBas[itk].Update();
768  mTEBas[itk].Invert();
769 
770  // Measured variables
771  for (int mem=0;mem<5;mem++) {
772  TString ts(TkPars::Name(mem)); ts+=itk;
773  AddPar(kMEAS,10*itk+mem,mTkDif[itk].Arr()+mem,1,ts,TkPars::Tiny(mem));
774  } }
775 
776  // Slack variables
777  AddPar(kSLAC,20,mLen,3,"Len" ,0.001);
778 
779  // Constrains
780  AddPar(kCNSR,30, mConr,7,"Same+Nrgy",1e-4);
781  Modify();
782  }
783  return TCFitData::Ready();
784 
785 }
786 //______________________________________________________________________________
787 void TCFitV0::Update()
788 {
789  if (!Modified()) return;
790 // Calc all constrains
791  double pos[3],dir[3];
792  TLorentzVector P4[3];
793  TVector3 Vx[3];
794  for (int itk=0;itk<2;itk++) {
795 // Calc all constrains
796  THelixTrack th;
797  mTkFit[itk] = mTkBas[itk];
798  mTkFit[itk]+= mTkDif[itk];
799  mTkFit[itk].Fill(th);
800  th.Eval(mLen[itk],pos,dir);
801  memcpy(mDConDL[itk],dir,sizeof(dir));
802  double P = mTkFit[itk].P();
803  double E = mTkFit[itk].E();
804  P4[itk].SetXYZT(dir[0]*P,dir[1]*P,dir[2]*P,E);
805  Vx[itk] = TVector3(pos);
806 // Calc DFcn
807  mTEBas[itk].Mpy(mTkDif[itk],mDFcn[itk]);
808 
809  }
810  P4[2]= P4[0]+P4[1];
811  (P4[2].Vect()*(-1.)).Unit().GetXYZ(mDConDL[2]);
812  Vx[2]= TVector3(mVx.x)+ P4[2].Vect().Unit()*mLen[2];
813  (Vx[0]-Vx[2]).GetXYZ(mConr+0);
814  (Vx[1]-Vx[2]).GetXYZ(mConr+3);
815  mConr[6] = P4[2].M() - mMas;
816  mTkFit[0].P4D(mP4d[0]);
817  mTkFit[1].P4D(mP4d[1]);
818 
819 
820  Modify(0);
821 }
822 //______________________________________________________________________________
823 int TCFitV0::Approx()
824 {
825 static int nCall=0; nCall++;
826 // Calc all constrains
827  double s[2],pos[3];
828  TVector3 Vx[3];
829  THelixTrack th[2];
830  for (int itk=0;itk<2;itk++) {mTkBas[itk].Fill(th[itk]);}
831 
832  double disMin=3e33;
833  for (int sgn=0;sgn<4;sgn++) {
834  if (sgn&1) th[0].Backward();
835  if (sgn&2) th[1].Backward();
836  s[0] = th[0].Path(th[1],&s[1]);
837  if (s[0]<100 && s[1]<100) {
838  for (int itk=0;itk<2;itk++) {
839  th[itk].Eval(s[itk],pos); Vx[itk].SetXYZ(pos[0],pos[1],pos[2]);}
840  double dis = (Vx[1]-Vx[0]).Mag();
841  if (dis < disMin) {
842  disMin=dis;
843  mLen[0] = (sgn&1)? -s[0]:s[0];
844  mLen[1] = (sgn&2)? -s[1]:s[1];
845  } }
846  if (sgn&1) th[0].Backward();
847  if (sgn&2) th[1].Backward();
848  }
849  if (disMin>100) return 1;
850  Vx[2]=(Vx[0]+Vx[1])*0.5;
851  mLen[2] = Vx[2].Mag();
852  return 0;
853 }
854 //______________________________________________________________________________
855 double TCFitV0::Fcn()
856 {
857  double chi2 = 0;
858 
859  for (int itk=0;itk<2;itk++) {
860 // Calc Fcn
861  chi2+=mTEBas[itk].Xi2(mTkDif[itk]);
862  }
863  return chi2;
864 }
865 //______________________________________________________________________________
866 double TCFitV0::Con(int ic)
867 {
868  return mConr[ic-30];
869 }
870 
871 //______________________________________________________________________________
872 double TCFitV0::DFcn(int ipar)
873 {
874  int itk = ipar/10;
875  assert( itk>=0 && itk <2);
876  int mem = ipar%10;
877  assert( mem>=0 && mem <5);
878  return mDFcn[itk][mem];
879 }
880 //______________________________________________________________________________
881 double TCFitV0::DDFcn(int ipar,int jpar)
882 {
883  int itk = ipar/10;
884  assert( itk>=0 && itk <2);
885  int jtk = jpar/10;
886  assert( jtk>=0 && jtk <2);
887  if (itk != jtk) return 0.;
888  int ivar = ipar%10;
889  int jvar = jpar%10;
890  assert(ivar<5 && jvar<5);
891  return mTEBas[itk].Get(ivar,jvar);
892 }
893 //______________________________________________________________________________
894 double TCFitV0::DCon(int icon,int ipar)
895 {
896  if (ipar <20) return TCFitData::DCon(icon,ipar);
897  if (icon==kCNRJ) return 0.;
898  int tcon = (icon-kCX_0)/3;
899  int jcon = (icon-kCX_0)%3;
900  int jpar = ipar-kLEN_0;
901  switch(10*tcon+jpar) {
902  case 00: ;case 11: ;
903  case 2: ;case 12: return mDConDL[jpar][jcon];
904  case 1: ;case 10: return 0.;
905  default: assert(0);
906  }
907  return 0;
908 }
909 //______________________________________________________________________________
910 void TCFitV0::Print(const char *name) const
911 {
912 static const char *ty[3]={"Meas","Slack","Cnsr"};
913 static const char *fx[2]={"Used","Fixed"};
914 
915  if (!name) name = "";
916  printf("TCFitV0::Print(%s) nMeas=%d(%d) nSlac=%d(%d) nCons=%d(%d)\n"
917  ,name
918  ,GetNMeas(),GetUMeas()
919  ,GetNSlac(),GetUSlac()
920  ,GetNCons(),GetUCons());
921  for (int i=0;i<kMaxId;i++) {
922  if (!fNams[i].Length()) continue;
923  const double *p = 0;
924  if (i< 5) { p = mTkFit[0].Arr()+i ;}
925  else if (i<15) { p = mTkFit[1].Arr()+i -10;}
926  if (!p) continue;
927  printf("%2d - %s\t",i,GetNam(i));
928  printf(" %s.%s ",ty[GetType(i)],fx[IsFixed(i)]);
929 
930  printf(" %g \n",*p);
931  }
932  TCFitData::Print(name);
933 }
934 //______________________________________________________________________________
935 #include "TRandom.h"
936 #include "TCanvas.h"
937 #include "TH1F.h"
938 #include "TSystem.h"
939 //______________________________________________________________________________
940 static double Pdk(double ma,double mb,double mc)
941 {
942  double ma2 = ma*ma,ma4=ma2*ma2;
943  double mb2 = mb*mb,mb4=mb2*mb2;
944  double mc2 = mc*mc,mc4=mc2*mc2;
945  double ans= (ma4+mb4+mc4-2*ma2*mb2-2*ma2*mc2-2*mb2*mc2);
946  if (ans<0.) ans = 0;
947  ans= sqrt(ans)/(2.*ma);
948  return ans;
949 }
950 //______________________________________________________________________________
951 void TCFitV0::Test(int mode)
952 {
953 static const double D0Mass=1.8646,PiMass=.13956995,KMass=.493677;
954 static const double D0Tau=0.415e-12,C=299792458e2;
955 static const double DX=0.003/*,DAng=0.003,DPti=0.1,DTan=0.005*/;
956 static const double D0P=1
957  ,D0Beta = D0P/sqrt(D0P*D0P+D0Mass*D0Mass)
958  ,D0Gama = 1./sqrt(1.-D0Beta*D0Beta)
959  ,D0Len = D0Tau*C*D0Gama*D0Beta;
960 static const double HZ= (0.000299792458 * 4.98478);
961 
962 static const int NHISTS = 4;
963 static TCanvas* myCanvas=0;
964 static TH1F *hh[]={0,0,0,0,0,0};
965 static const char *hNams[]= {"dL","dL/L", "Xi2","Mass-D0"};
966 static const double LOW[] = {-D0Len,-3, 0.,-1.5};
967 static const double UPP[] = {-D0Len, 3,10., 1.5};
968 
969 
970 
971  if(!myCanvas) myCanvas=new TCanvas("TCircleFitter_Test","",600,800);
972  myCanvas->Clear();
973  myCanvas->Divide(1,NHISTS);
974  for (int i=0;i<NHISTS;i++) {
975  delete hh[i]; hh[i]= new TH1F(hNams[i],hNams[i],100,LOW[i],UPP[i]);
976  myCanvas->cd(i+1); hh[i]->Draw();
977  }
978 
979 TLorentzVector p4[3];
980 TVector3 dcayPos,d0Mom;
981 THelixTrack Khlx,Pihlx;
982 double Q = Pdk(D0Mass,PiMass,KMass);
983 
984  TCFitV0 dat;
985  TCFit tc("TestD0",&dat);
986  dat.mTkBas[0].SetHz(1.);
987  dat.mTkBas[1].SetHz(1.);
988 
989  for (int iev=0; iev <500; iev++) {
990  tc.Reset(); dat.Reset();
991  double dist = gRandom->Exp(D0Len);
992  double wk[9];
993  gRandom->Sphere(wk[0],wk[1],wk[2],dist);
994  dcayPos.SetXYZ(wk[0],wk[1],wk[2]);
995  d0Mom = dcayPos.Unit()*D0P;
996  p4[2].SetVectM( d0Mom,D0Mass);
997  gRandom->Sphere(wk[0],wk[1],wk[2],Q);
998  p4[0].SetXYZM(-wk[0],-wk[1],-wk[2],KMass );
999  p4[1].SetXYZM( wk[0], wk[1], wk[2],PiMass);
1000  TVector3 boost = p4[2].BoostVector();
1001  dat.mMas = D0Mass;
1002  for (int ip=0;ip<2;ip++) {
1003  p4[ip].Boost(boost);
1004  p4[ip].GetXYZT(wk+3);
1005  dcayPos.GetXYZ (wk+0);
1006  double ptin = 1./p4[ip].Pt();
1007  if (ip) ptin=-ptin;
1008 
1009  THelixTrack th(wk+0,wk+3,ptin*HZ);
1010  th.Backward();
1011  double s =th.Path(0.,0.);
1012  th.Move(s);th.Backward();
1013  th.Eval(0.,wk,wk+3);
1014  TVector3 pos(wk),dir(wk+3);
1015  dat.mTkBas[ip].Reset();
1016  dat.mTkBas[ip].Set(pos,dir,ptin);
1017  dat.mTkBas[ip].mass= (ip)? PiMass:KMass;
1018  dat.mTkBas[ip].Update();
1019  dat.mTEBas[ip].Set(0,0,pow(DX,2));
1020  dat.mTEBas[ip].Set(1,1,pow(DX,2));
1021  dat.mTkBas[ip].Rand(dat.mTEBas[ip]);
1022  }
1023  dat.Ready();
1024  if (mode==1) dat.FixPar(kCNRJ);
1025  if (tc.Fit()) continue;
1026  double dL = dat.GetPar(kLEN_2)-dist;
1027  double eL = sqrt(dat.ErMx(kLEN_2,kLEN_2));
1028  hh[0]->Fill(dL/eL);
1029  hh[1]->Fill(dL/dist);
1030  hh[2]->Fill(dat.GetFcn()/dat.GetNDF());
1031  TLorentzVector D0V = dat.mTkFit[0].P4()+dat.mTkFit[1].P4();
1032  hh[3]->Fill(D0V.M()-D0Mass);
1033 
1034  }
1035  myCanvas->Modified();
1036  myCanvas->Update();
1037  // while(!gSystem->ProcessEvents()){};
1038 }
1039 
THelixTrack::Eval
double Eval(double step, double *xyz, double *dir, double &rho) const
Evaluate params with given step along helix.
Definition: THelixTrack.cxx:1063
TCFitData::fModi
int fModi
fail flag, fit is impossible
Definition: TCFit.h:169
TkPars::tanl
double tanl
tangent of the track momentum dip angle
Definition: TCFit.h:228
TCFitData::fNPars
int fNPars[3]
&amp;1=1 fcn calculated from error matrix
Definition: TCFit.h:171
TkPars::mass
double mass
Mass of track.
Definition: TCFit.h:234
TkPars
Definition: TCFit.h:196
TCFit
Definition: TCFit.h:19
TCFitData::fTyps
short fTyps[kMaxId+1]
Tiny values still modifying fcn.
Definition: TCFit.h:176
TCFitV0
Definition: TCFit.h:267
TkPars::curv
double curv
signed curvature [sign = sign(-qB)]
Definition: TCFit.h:230
TCL::trsa
static float * trsa(const float *s, const float *a, float *b, int m, int n)
Definition: TCernLib.cxx:1111
TkPars::hz
double hz
Z component magnetic field in units Pt(Gev) = Hz * RCurv(cm)
Definition: TCFit.h:232
THelixTrack::Backward
void Backward()
Change direction.
Definition: THelixTrack.cxx:369
TkErrs
Definition: TCFit.h:237
TkPars::dca
double dca
point
Definition: TCFit.h:222
TCL::trsinv
static float * trsinv(const float *g, float *gi, int n)
Definition: TCernLib.cxx:1160
TCL::trasat
static float * trasat(const float *a, const float *s, float *r, int m, int n)
Definition: TCernLib.cxx:540
TkPars::phi
double phi
angle between track direction and X-axis in xy plane
Definition: TCFit.h:224
TCFitData::fPars
double * fPars[kMaxId+1]
number of &quot;slack&quot; parameters
Definition: TCFit.h:174
TCFitData::fEnd
char fEnd[1]
Current value of fcn [1].
Definition: TCFit.h:181
TkPars::ptin
double ptin
signed invert pt [sign = sign(-qB)]
Definition: TCFit.h:226
TCFitData::fNFixs
int fNFixs[3]
number of &quot;measured&quot; ,slack, constrains
Definition: TCFit.h:172
TEST0
Definition: TCFit.cxx:245
THelixTrack::Move
double Move(double step)
Move along helix.
Definition: THelixTrack.cxx:452
TCFitData::AddPar
int AddPar(int tyPar, int idPar, double *par, int nPars=1, const char *name="", double tiny=0.)
Definition: TCFit.cxx:416
Generated by   doxygen 1.8.5