00001
00002
00003
00004
00005 #include <stdio.h>
00006 #include <stdlib.h>
00007 #include <assert.h>
00008 #include <math.h>
00009 #include "TError.h"
00010 #include "TCFit.h"
00011 #include "TNumDeriv.h"
00012 #include "TVector3.h"
00013 #include "TRandom.h"
00014 #include "TLorentzVector.h"
00015 #include "THelixTrack.h"
00016 #if ROOT_VERSION_CODE < 331013
00017 #include "TCL.h"
00018 #else
00019 #include "TCernLib.h"
00020 #endif
00021 ClassImp(TCFit)
00022
00023 TCFit::TCFit(const char *name,TCFitData *dat) :TNamed(name,"")
00024 {
00025 fDebug = 1;
00026 SetData(dat);
00027 for (int i=0;i<5;i++) (&fBigM)[i]=new TMatrixD;
00028
00029 }
00030
00031 void TCFit::Reset()
00032 {
00033 for (int i=0;i<5;i++) {*((&fBigM)[i])=0.;}
00034 fIter =0;
00035 fCuts =0;
00036 fMaxIter =500;
00037 fMaxCuts =6;
00038 fUPars = 0;
00039 fUMeas = 0;
00040 fUCons = 0;
00041 }
00042
00043 TCFit::~TCFit()
00044 {
00045 for (int i=0;i<4;i++) delete (&fBigM)[i];
00046 }
00047
00048 int TCFit::SetData(TCFitData *dat)
00049 {
00050 fDat = dat;
00051 fDat->SetFitter(this);
00052 return 0;
00053 }
00054
00055 int TCFit::CheckOut()
00056 {
00057 int nBig = fUPars+fUCons;
00058
00059 fDat->Evaluate();
00060
00061 if (!fIter) {
00062 fBigM->ResizeTo(nBig,nBig); (*fBigM)*=0.0;
00063 fBigI->ResizeTo(nBig,nBig);
00064 fBigB->ResizeTo(nBig,1);
00065 fOldP->ResizeTo(nBig,1);
00066 fAddP->ResizeTo(nBig,1);
00067 }
00068 for (int j1=0;j1<fUMeas; j1++) {
00069 int i1 = fDat->GetId(j1);
00070 if (!fIter) (*fOldP )[j1][ 0] = fDat->GetPar(i1);
00071 (*fBigB )[j1][ 0] = -fDat->DFcn(i1);
00072 (*fBigM)[j1][j1] = fDat->DDFcn(i1,i1);
00073 for (int j2=0;j2<j1; j2++) {
00074 int i2 = fDat->GetId(j2);
00075 (*fBigM)[j1][j2] = fDat->DDFcn(i1,i2);
00076 (*fBigM)[j2][j1] = (*fBigM)[j1][j2];
00077 }
00078 }
00079 for (int jc=0;jc<fUCons; jc++) {
00080 int ic = fDat->GetId(jc+fUPars);
00081 double cc = fDat->GetPar(ic);
00082 (*fBigB )[jc+fUPars][0] = -cc;
00083 for (int jx=0;jx<fUPars; jx++) {
00084 int ix = fDat->GetId(jx);
00085 double cx = fDat->DCon(ic,ix);
00086 (*fBigM)[fUPars+jc][jx] = cx;
00087 (*fBigM)[jx][fUPars+jc] = cx;
00088 } }
00089
00090
00091 for (int jx=0; jx<fUPars;jx++) {
00092 for (int jc=0;jc<fUCons;jc++) {
00093 (*fBigB )[jx][ 0] -= (*fOldP)[fUPars+jc][0]*(*fBigM)[jx][fUPars+jc];
00094 } }
00095
00096 if (fDebug>=3) {
00097 fBigM->Print("CheckOut.BigM");
00098 fBigB->Print("CheckOut.BigB");
00099 fDat->Print("CheckOut.dat");
00100 }
00101
00102
00103 return 0;
00104 }
00105
00106 int TCFit::CheckIn()
00107 {
00108 for (int jx=0;jx<fUPars;jx++) {
00109 int ix = fDat->GetId(jx);
00110 fDat->GetPar(ix) = (*fOldP)[jx][0] + (*fAddP)[jx][0];
00111 }
00112 fDat->Evaluate();
00113
00114 return 0;
00115 }
00116
00117 int TCFit::Fit()
00118 {
00119 if (fDat->Ready()) return kBadAprx;
00120 fUPars = fDat->GetUPars();
00121 fUMeas = fDat->GetUMeas();
00122 fUCons = fDat->GetUCons();
00123
00124 fCuts = 0; fFitRes=0;
00125
00126 fFcnQA[0]=3e33; fFcnQA[1]=3e33;
00127 fConQA[0]=3e33; fConQA[1]=3e33;
00128 fAddQA[0]=3e33; fAddQA[1]=3e33;
00129
00130
00131
00132
00133 for (fIter=0; fIter < fMaxIter; fIter++)
00134 {
00135 fAkt = CheckStep();
00136 switch (fAkt) {
00137
00138 case kNextStep: (*fOldP)+=(*fAddP); FitStep(); break;
00139
00140 case kNextCut: CutStep(); break;
00141
00142 case kEndFit: EndStep(); return fFitRes;
00143
00144 case kBadFit: EndStep(); return fFitRes;
00145
00146 default: assert(0);
00147 }
00148 if (fDebug>=2) PriStep("OneIter");
00149
00150 }
00151 fFitRes=kTooIter;
00152 EndStep();
00153 return fFitRes;
00154 }
00155
00156 int TCFit::FitStep()
00157 {
00158 CheckOut();
00159 fFcnQA[1]=fFcnQA[0];
00160 fConQA[1]=fConQA[0];
00161 fAddQA[1]=fAddQA[0];
00162
00163 (*fBigI) = (*fBigM);
00164 double det=0;
00165 (*fBigI).Invert(&det);
00166 (*fAddP) = (*fBigI)*(*fBigB);
00167
00168 return 0;
00169 }
00170
00171 int TCFit::CutStep()
00172 {
00173 fCuts++; fIter--;
00174 (*fAddP)*=0.5; return 0;
00175 }
00176
00177 int TCFit::EndStep()
00178 {
00179 fDat->SetFail(fFitRes);
00180 if (fDebug>=1) PriStep("EndStep");
00181 if (fDebug>=2) fDat->Print("EndStep");
00182 return 0;
00183 }
00184
00185 int TCFit::CheckStep()
00186 {
00187 fFitRes = 0;
00188
00189 if (fIter==0) return kNextStep;
00190 CheckIn();
00191 fFcnQA[0] = fDat->GetFcn();
00192 fConQA[0] = 0;
00193 fAddQA[0] = 0;
00194 for (int jc=0;jc<fUCons; jc++) {
00195 int ic = fDat->GetId(jc+fUPars);
00196 fConQA[0]+= fabs(fDat->GetPar(ic)/fDat->GetTiny(ic));
00197 }
00198 if (fUCons>1) fConQA[0]/=fUCons;
00199 for (int jx=0;jx<fUMeas; jx++) {
00200 int ix = fDat->GetId(jx);
00201 fAddQA[0]+= fabs((*fAddP)[jx][0]/fDat->GetTiny(ix));
00202 }
00203 fAddQA[0]/=fUMeas;
00204
00205 if (fAddQA[0]<1. && fConQA[0] <1.) return kEndFit;
00206
00207 if (fFcnQA[0] > fDat->GetBigFcn() ) fFitRes |=kBadFcn;
00208 if (fConQA[0] > 1e10 ) fFitRes |=kBadCon;
00209 if (fFitRes) return kBadFit;
00210
00211 if (fAddQA[0]<fAddQA[1]
00212 || fConQA[0]<fConQA[1]
00213 || fFcnQA[0]<fFcnQA[1]
00214 || fCuts>fMaxCuts) {
00215 fCuts = 0;
00216 fAddQA[1]=fAddQA[0];
00217 fConQA[1]=fConQA[0];
00218 fFcnQA[1]=fFcnQA[0];
00219 return kNextStep;
00220 }
00221 return kNextCut;
00222 }
00223
00224 double TCFit::ErMx(int jcol,int jrow) const
00225 {
00226 return (*fBigI)[jcol][jrow];
00227 }
00228
00229 int TCFit::PriStep(const char *where)
00230 {
00231 if (where==0)where="";
00232 printf("PriStep(%s) Iter=%d Cut=%d Akt=%d Fcn=%g(%g) Con=%g(%g) Add=%g(%g)\n"
00233 ,where,fIter,fCuts,fAkt
00234 ,fFcnQA[0],fFcnQA[1],fConQA[0],fConQA[1],fAddQA[0],fAddQA[1]);
00235
00236 return 0;
00237 }
00238
00239
00240
00241
00242
00243
00244
00245 class TEST0 : public TCFitData {
00246 public:
00247 TEST0();
00248
00249 double Fcn();
00250 double Con(int icon);
00251
00252 double MyDFcn(int ix);
00253 double MyDDFcn(int ix);
00254 void Update() {}
00255 int Approx() {return 0;}
00256 public:
00257 double mX[3],mC;
00258
00259 };
00260
00261 TEST0::TEST0():TCFitData("TEST0")
00262 {
00263 mX[0]=mX[1]=mX[2]=1.1;
00264 mC=0;
00265 AddPar(0,0, mX,3,"X",0.0001);
00266 AddPar(2,3,&mC,1,"C",0.0001);
00267 }
00268
00269 double TEST0::Fcn()
00270 {
00271 return mX[0]*mX[0]/1. + mX[1]*mX[1]/4 + mX[2]*mX[2]/9;
00272 }
00273
00274 double TEST0::MyDFcn(int ix)
00275 {
00276 switch (ix) {
00277 case 0: return 2*mX[0];
00278 case 1: return 2*mX[1]/4;
00279 case 2: return 2*mX[2]/9;
00280 }
00281 assert(0); return 0.;
00282
00283 }
00284
00285 double TEST0::MyDDFcn(int ix)
00286 {
00287 switch (ix) {
00288 case 0: return 2.;
00289 case 1: return 2./4;
00290 case 2: return 2./9;
00291 }
00292 assert(0); return 0.;
00293
00294 }
00295
00296 double TEST0::Con(int)
00297 {
00298 return mX[2]-1.;
00299 }
00300
00301
00302
00303 void TCFit::Test0()
00304 {
00305 TEST0 t0;
00306 TCFit f0("Test0",&t0);
00307 for (int itst=0;itst<3;itst++) {
00308 double d1 = t0.DFcn(itst);
00309 double d2 = t0.MyDFcn(itst);
00310 double dd1 = t0.DDFcn(itst,itst);
00311 double dd2 = t0.MyDDFcn(itst);
00312 printf("%d %g %g %g %g\n",itst,d1,d2,dd1,dd2);
00313 }
00314 f0.Fit();
00315 }
00316
00317
00318
00319 ClassImp(TCFitData)
00320
00321 class Deriv1st : public TNumDeriv {
00322 public:
00323 Deriv1st(TCFitData* dat);
00324 Double_t Fcn(Double_t arg);
00325 void SetKind(int kind) { fKind=kind;}
00326 private:
00327 TCFitData *fFitData;
00328 int fKind;
00329 };
00330
00331 Deriv1st::Deriv1st(TCFitData* dat):TNumDeriv("Deriv1st")
00332 {
00333 fFitData=dat;
00334 fKind = -1;
00335 }
00336
00337 Double_t Deriv1st::Fcn(Double_t arg)
00338 {
00339 double save = fFitData->GetPar(fIArg);
00340 SetStep(fFitData->GetTiny(fIArg));
00341 fFitData->GetPar(fIArg)+=arg;
00342 fFitData->Update(); fFitData->Modify(0);
00343 double ans = (fKind==(-1)) ? fFitData->Fcn() :fFitData->Con(fKind);
00344 fFitData->GetPar(fIArg)=save;
00345 fFitData->Update(); fFitData->Modify(0);
00346 return ans;
00347 }
00348
00349
00350
00351 class Deriv2nd : public TNumDeriv {
00352 public:
00353 Deriv2nd(TCFitData* dat);
00354 ~Deriv2nd();
00355 Double_t Fcn(Double_t arg);
00356 void SetKind(int kind) { fD1->SetKind(kind);}
00357 void SetIJArg(int i,int j);
00358 private:
00359 TCFitData *fFitData;
00360 Deriv1st *fD1;
00361 };
00362
00363 Deriv2nd::Deriv2nd(TCFitData* dat):TNumDeriv("Deriv2nd")
00364 {
00365 fFitData = dat;
00366 fD1 = new Deriv1st(fFitData);
00367 }
00368
00369 Deriv2nd::~Deriv2nd()
00370 {
00371 delete fD1;
00372 }
00373
00374 void Deriv2nd::SetIJArg(int i,int j)
00375 {
00376 this->SetIArg(i);
00377 fD1 ->SetIArg(j);
00378 }
00379
00380 Double_t Deriv2nd::Fcn(Double_t arg)
00381 {
00382 double save = fFitData->GetPar(fIArg);
00383 fFitData->GetPar(fIArg)+=arg;
00384 fFitData->Update(); fFitData->Modify(0);
00385 double ans = fD1->DFcn(0.);
00386 fFitData->GetPar(fIArg)=save;
00387 fFitData->Update(); fFitData->Modify(0);
00388 return ans;
00389 }
00390
00391
00392 TCFitData::TCFitData(const char *name, const char *title):TNamed(name,title)
00393 {
00394
00395 Reset();
00396 fD1st = new Deriv1st(this);
00397 fD2nd = new Deriv2nd(this);
00398 }
00399
00400 void TCFitData::Reset()
00401 {
00402 memset(fBeg,0,fEnd-fBeg+1);
00403 memset(fFixs,1,sizeof(fFixs));
00404 for (int i=0; i<kMaxId;i++) fNams[i]="";
00405 fFcn[1]=1e-6;
00406 fFcn[2]=1e+6;
00407 }
00408
00409
00410 TCFitData::~TCFitData()
00411 {
00412 delete fD1st;
00413 delete fD2nd;
00414 }
00415
00416 int TCFitData::AddPar(int tyPar, int idPar,double *par,int nPars,const char *name,double tiny)
00417 {
00418 assert(0<=tyPar && tyPar<=2);
00419 assert(0<=idPar && idPar+nPars<kMaxId);
00420 if (fFail ) return fFail;
00421 if (!name) name="";
00422 if (tiny<=0) tiny = 1e-3;
00423
00424 for (int j=0;j<nPars;j++) {
00425 TString ts;
00426 if (name[0]) ts=name; else ts=idPar;
00427 if (nPars>1) { ts +="["; ts +=j; ts+="]";}
00428 fNams[idPar+j]= ts;
00429 fTiny[idPar+j]= tiny;
00430 fPars[idPar+j]= par+j;
00431 fFixs[idPar+j]= 0;
00432 fTyps[idPar+j]= tyPar;
00433 }
00434 fNPars[tyPar]+=nPars;
00435 assert(fNPars[tyPar]<=kMaxId);
00436 return 0;
00437 }
00438
00439 Int_t TCFitData::GetId(const char *name) const
00440 {
00441 int i=0;
00442 int n = fNPars[0]+fNPars[1]+fNPars[2];
00443 for (int j=0;i<n;j++) {
00444 if (!fNams[j].Length()) continue;
00445 i++;
00446 if (fNams[j]==name) return j;
00447 }
00448 printf("TCFitData::GetId(\"%s\") UNKNOWN name\n",name);
00449 return -1;
00450 }
00451
00452 const char *TCFitData::GetNam(Int_t idx) const
00453 {
00454 return fNams[idx].Data();
00455 }
00456
00457 int TCFitData::GetType(int id) const
00458 {
00459
00460 return fTyps[id];
00461 }
00462
00463
00464
00465 int TCFitData::GetId(int jdx) const
00466 {
00467 return fIndx[jdx];
00468 }
00469
00470 int TCFitData::GetJd(int idx) const
00471 {
00472 return fJndx[idx];
00473 }
00474
00475
00476 void TCFitData::FixPar(int idx,int yes)
00477 {
00478 yes = (yes!=0);
00479 assert(fNams[idx].Length());
00480 if (fFixs[idx]==yes) return;
00481 fFixs[idx]=yes;
00482 fNFixs[fTyps[idx]]+= ( (yes)? 1:-1 );
00483
00484 }
00485
00486 int TCFitData::IsFixed(int idx) const
00487 { return fFixs[idx];}
00488
00489
00490 double TCFitData::GetPar (int idx) const
00491 {return *fPars[idx];}
00492
00493 double &TCFitData::GetPar (int idx)
00494 {
00495 fModi=1; return *fPars[idx];
00496 }
00497
00498 void TCFitData::Evaluate()
00499 {
00500 if (Modified()) Update(); Modify(0);
00501 fFcn[0] = Fcn();
00502 for (int jc=GetUPars();jc<GetUPars()+GetUCons();jc++) {
00503 int ic = GetId(jc);
00504 GetPar(ic) = Con(ic);
00505 }
00506 }
00507
00508 int TCFitData::Ready()
00509 {
00510 int n = 0;
00511 for (int ity=0;ity<3;ity++) {
00512 for (int i=0;i<kMaxId;i++) {
00513 if (!fNams[i].Length()) continue;
00514 if ( fFixs[i]) continue;
00515 if ( fTyps[i]!=ity) continue;
00516 fIndx[n]=i;
00517 fJndx[i]=n++;
00518 } }
00519 assert(n==fNPars[0]-fNFixs[0]+fNPars[1]-fNFixs[1]+fNPars[2]-fNFixs[2]);
00520 return Approx();
00521 }
00522
00523
00524 Double_t TCFitData::DFcn(int ipar)
00525 {
00526 fD1st->SetKind(-1);
00527 fD1st->SetIArg(ipar);
00528 return fD1st->DFcn();
00529 }
00530
00531
00532 Double_t TCFitData::DDFcn(int ipar,int jpar)
00533 {
00534 fD2nd->SetKind(-1);
00535 fD2nd->SetIJArg(ipar,jpar);
00536 return fD2nd->DFcn();
00537 }
00538
00539
00540 Double_t TCFitData::Con(int)
00541 {
00542 assert(0);
00543 }
00544
00545
00546 Double_t TCFitData::DCon(int icon,int ipar)
00547 {
00548 fD1st->SetKind(icon);
00549 fD1st->SetIArg(ipar);
00550 return fD1st->DFcn();
00551 }
00552
00553 double TCFitData::ErMx(int icol,int irow) const
00554 {
00555 int jcol = GetJd(icol);
00556 int jrow = GetJd(irow);
00557 return fFitter->ErMx(jcol,jrow);
00558 }
00559
00560 void TCFitData::Print(const char *name) const
00561 {
00562 static const char *ty[3]={"Meas","Slack","Cnsr"};
00563 static const char *fx[2]={"Used","Fixed"};
00564
00565 if (!name) name = "";
00566 printf("TCFitData::Print(%s) nMeas=%d(%d) nSlac=%d(%d) nCons=%d(%d)\n"
00567 ,name
00568 ,GetNMeas(),GetUMeas()
00569 ,GetNSlac(),GetUSlac()
00570 ,GetNCons(),GetUCons());
00571 for (int i=0;i<kMaxId;i++) {
00572 if (!fNams[i].Length()) continue;
00573 printf("%2d - %s\t",i,fNams[i].Data());
00574 printf(" %s.%s ",ty[fTyps[i]],fx[fFixs[i]]);
00575 printf(" %g \n",*fPars[i]);
00576 }
00577
00578 }
00579
00580
00581
00582 void TkPars::Print(const char *tit) const
00583 {
00584 if (!tit) tit="";
00585 printf("TkPars::Print(%s) \tDca=%g Z=%g Phi=%g ptin=%g tanl=%g curv=%g"
00586 ,tit,dca,z,phi,ptin,tanl,curv);
00587
00588 }
00589
00590 void TkPars::Reset()
00591 {
00592 double tmp = hz;
00593 memset(this,0,sizeof(TkPars));
00594 hz = tmp;
00595 }
00596
00597 void TkPars::SetHz(double fact)
00598 {
00599 hz = (0.000299792458 * 4.98478)*fact;
00600 }
00601
00602 TkPars &TkPars::operator+=(const TkPars &a)
00603 {
00604 for (int i=0;i<5;i++) Arr()[i]+=a.Arr()[i];
00605 if (phi<-M_PI) phi+=M_PI*2;
00606 if (phi> M_PI) phi-=M_PI*2;
00607 return *this;
00608 }
00609
00610 TLorentzVector TkPars::P4() const
00611 {
00612 double pt = fabs(1./ptin);
00613 double e = sqrt(pt*pt*(1.+tanl*tanl)+mass*mass);
00614 TLorentzVector v(pt*cos(phi),pt*sin(phi),pt*tanl,e);
00615 return v;
00616 }
00617
00618 void TkPars::P4D(double D[4][5]) const
00619 {
00620 enum {Kdca,Kz,Kphi,Kptin,Ktanl};
00621
00622 memset(D[0],0,4*5*sizeof(D[0][0]));
00623 double pt = fabs(1./ptin);
00624 double dpt = -pt/ptin;
00625 double e = sqrt(pt*pt*(1.+tanl*tanl)+mass*mass);
00626
00627
00628 D[3][Kptin] = pt*(1.+tanl*tanl)/e*dpt;
00629 D[3][Ktanl] = pt*pt*tanl/e*dpt;
00630
00631 D[0][Kptin] = cos(phi)*dpt;
00632 D[1][Kptin] = sin(phi)*dpt;
00633 D[2][Kptin] = tanl*dpt;
00634 D[2][Ktanl] = pt;
00635
00636 }
00637
00638 void TkPars::Set(const TVector3 &v3,const TVector3 &d3,double pti)
00639 {
00640 dca = v3.Perp();
00641 phi = d3.Phi();
00642 z = v3[2];
00643 if ((v3.Cross(d3)).Z()<0) dca=-dca;
00644 ptin = pti;
00645 tanl = d3.Pz()/d3.Pt();
00646 }
00647
00648
00649 void TkPars::Get(TVector3 *v3,TVector3 *d3,double *pts) const
00650 {
00651 if (v3 ) { v3->SetXYZ(dca*sin(phi),-dca*cos(phi),z) ;}
00652 if (d3 ) { v3->SetXYZ(cos(phi),sin(phi),tanl); v3->SetMag(1.);}
00653 if (pts) {*pts = 1./ptin ;}
00654 }
00655
00656 TVector3 TkPars::V3() const
00657 {
00658 TVector3 v3; Get(&v3); return v3;
00659 }
00660
00661 void TkPars::Fill(THelixTrack &hlx)
00662 {
00663 TLorentzVector p4 = P4();
00664 TVector3 v3 = V3();
00665 double h[6];
00666 v3.GetXYZ(h);
00667 p4.Vect().GetXYZ(h+3);
00668 hlx.Set(h,h+3,curv);
00669
00670 }
00671
00672 const char* TkPars::Name(int mem)
00673 {
00674 static const char* nams[]={ "Dca", "Z", "Phi", "Ptin", "TanL"};
00675 return nams[mem];
00676 }
00677
00678 double TkPars::Tiny(int mem)
00679 {
00680 static const double tiny[]={3.14/180/10., 0.01, 3.14/180/10, 0.01, 0.001};
00681 return tiny[mem];
00682 }
00683
00684 void TkPars::Rand(const TkErrs &errs)
00685 {
00686 for (int iv=0;iv<5;iv++) {
00687 Arr()[iv]+= gRandom->Gaus(0.,sqrt(errs.Get(iv,iv)));
00688 }
00689 Update();
00690 }
00691
00692
00693
00694 void TkErrs::Reset()
00695 {
00696 static const double DX=0.3,DZ=0.3,DPhi=0.03,DPti=1.,DTan=0.05;
00697 memset(emx,0,sizeof(emx));
00698 double fak = 0.3*0.3;
00699 emx[ 0]=DX*DX *fak;
00700 emx[ 2]=DZ*DZ *fak;
00701 emx[ 5]=DPhi*DPhi *fak;
00702 emx[ 9]=DPti*DPti *fak;
00703 emx[14]=DTan*DTan *fak;
00704 }
00705
00706 void TkErrs::Set(int i,int j,double err)
00707 {
00708 if (i<j) { int ii=i; i=j; j=ii;}
00709 emx[((i+1)*i)/2+j]=err;
00710 }
00711
00712 double TkErrs::Get(int i,int j) const
00713 {
00714 if (i<j) { int ii=i; i=j; j=ii;}
00715 return emx[((i+1)*i)/2+j];
00716 }
00717
00718 void TkErrs::Invert()
00719 {
00720 TCL::trsinv(emx,emx,5);
00721 }
00722
00723 double TkErrs::Xi2(const TkPars &pars) const
00724 {
00725 double chi2;
00726 TCL::trasat(pars.Arr(),emx,&chi2,1,5);
00727 return chi2;
00728 }
00729
00730 void TkErrs::Mpy(const TkPars &pars,double der[5]) const
00731 {
00732 TCL::trsa(emx, pars.Arr() ,der,5,1);
00733 }
00734
00735
00736 TVector3 VxPars::V3() const
00737 { return TVector3(x); }
00738
00739
00740
00741 ClassImp(TCFitV0)
00742
00743 TCFitV0::TCFitV0():TCFitData("TCFitV0","Two prong decay fit")
00744 {
00745 Reset();
00746 }
00747
00748 void TCFitV0::Reset()
00749 {
00750 for (int i=0;i<2;i++) {
00751 mTkBas[i].Reset();
00752 mTEBas[i].Reset();
00753 mTkFit[i].Reset();
00754 mTkDif[i].Reset();
00755 }
00756 memset(mBeg,0,mEnd-mBeg+1);
00757 TCFitData::Reset();
00758 }
00759
00760 int TCFitV0::Ready()
00761 {
00762 if (!mReady) {
00763 mReady=46;
00764 TCFitData::Reset();
00765
00766 for (int itk=0;itk<2;itk++) {
00767 mTkBas[itk].Update();
00768 mTEBas[itk].Invert();
00769
00770
00771 for (int mem=0;mem<5;mem++) {
00772 TString ts(TkPars::Name(mem)); ts+=itk;
00773 AddPar(kMEAS,10*itk+mem,mTkDif[itk].Arr()+mem,1,ts,TkPars::Tiny(mem));
00774 } }
00775
00776
00777 AddPar(kSLAC,20,mLen,3,"Len" ,0.001);
00778
00779
00780 AddPar(kCNSR,30, mConr,7,"Same+Nrgy",1e-4);
00781 Modify();
00782 }
00783 return TCFitData::Ready();
00784
00785 }
00786
00787 void TCFitV0::Update()
00788 {
00789 if (!Modified()) return;
00790
00791 double pos[3],dir[3];
00792 TLorentzVector P4[3];
00793 TVector3 Vx[3];
00794 for (int itk=0;itk<2;itk++) {
00795
00796 THelixTrack th;
00797 mTkFit[itk] = mTkBas[itk];
00798 mTkFit[itk]+= mTkDif[itk];
00799 mTkFit[itk].Fill(th);
00800 th.Eval(mLen[itk],pos,dir);
00801 memcpy(mDConDL[itk],dir,sizeof(dir));
00802 double P = mTkFit[itk].P();
00803 double E = mTkFit[itk].E();
00804 P4[itk].SetXYZT(dir[0]*P,dir[1]*P,dir[2]*P,E);
00805 Vx[itk] = TVector3(pos);
00806
00807 mTEBas[itk].Mpy(mTkDif[itk],mDFcn[itk]);
00808
00809 }
00810 P4[2]= P4[0]+P4[1];
00811 (P4[2].Vect()*(-1.)).Unit().GetXYZ(mDConDL[2]);
00812 Vx[2]= TVector3(mVx.x)+ P4[2].Vect().Unit()*mLen[2];
00813 (Vx[0]-Vx[2]).GetXYZ(mConr+0);
00814 (Vx[1]-Vx[2]).GetXYZ(mConr+3);
00815 mConr[6] = P4[2].M() - mMas;
00816 mTkFit[0].P4D(mP4d[0]);
00817 mTkFit[1].P4D(mP4d[1]);
00818
00819
00820 Modify(0);
00821 }
00822
00823 int TCFitV0::Approx()
00824 {
00825 static int nCall=0; nCall++;
00826
00827 double s[2],pos[3];
00828 TVector3 Vx[3];
00829 THelixTrack th[2];
00830 for (int itk=0;itk<2;itk++) {mTkBas[itk].Fill(th[itk]);}
00831
00832 double disMin=3e33;
00833 for (int sgn=0;sgn<4;sgn++) {
00834 if (sgn&1) th[0].Backward();
00835 if (sgn&2) th[1].Backward();
00836 s[0] = th[0].Path(th[1],&s[1]);
00837 if (s[0]<100 && s[1]<100) {
00838 for (int itk=0;itk<2;itk++) {
00839 th[itk].Eval(s[itk],pos); Vx[itk].SetXYZ(pos[0],pos[1],pos[2]);}
00840 double dis = (Vx[1]-Vx[0]).Mag();
00841 if (dis < disMin) {
00842 disMin=dis;
00843 mLen[0] = (sgn&1)? -s[0]:s[0];
00844 mLen[1] = (sgn&2)? -s[1]:s[1];
00845 } }
00846 if (sgn&1) th[0].Backward();
00847 if (sgn&2) th[1].Backward();
00848 }
00849 if (disMin>100) return 1;
00850 Vx[2]=(Vx[0]+Vx[1])*0.5;
00851 mLen[2] = Vx[2].Mag();
00852 return 0;
00853 }
00854
00855 double TCFitV0::Fcn()
00856 {
00857 double chi2 = 0;
00858
00859 for (int itk=0;itk<2;itk++) {
00860
00861 chi2+=mTEBas[itk].Xi2(mTkDif[itk]);
00862 }
00863 return chi2;
00864 }
00865
00866 double TCFitV0::Con(int ic)
00867 {
00868 return mConr[ic-30];
00869 }
00870
00871
00872 double TCFitV0::DFcn(int ipar)
00873 {
00874 int itk = ipar/10;
00875 assert( itk>=0 && itk <2);
00876 int mem = ipar%10;
00877 assert( mem>=0 && mem <5);
00878 return mDFcn[itk][mem];
00879 }
00880
00881 double TCFitV0::DDFcn(int ipar,int jpar)
00882 {
00883 int itk = ipar/10;
00884 assert( itk>=0 && itk <2);
00885 int jtk = jpar/10;
00886 assert( jtk>=0 && jtk <2);
00887 if (itk != jtk) return 0.;
00888 int ivar = ipar%10;
00889 int jvar = jpar%10;
00890 assert(ivar<5 && jvar<5);
00891 return mTEBas[itk].Get(ivar,jvar);
00892 }
00893
00894 double TCFitV0::DCon(int icon,int ipar)
00895 {
00896 if (ipar <20) return TCFitData::DCon(icon,ipar);
00897 if (icon==kCNRJ) return 0.;
00898 int tcon = (icon-kCX_0)/3;
00899 int jcon = (icon-kCX_0)%3;
00900 int jpar = ipar-kLEN_0;
00901 switch(10*tcon+jpar) {
00902 case 00: ;case 11: ;
00903 case 2: ;case 12: return mDConDL[jpar][jcon];
00904 case 1: ;case 10: return 0.;
00905 default: assert(0);
00906 }
00907 return 0;
00908 }
00909
00910 void TCFitV0::Print(const char *name) const
00911 {
00912 static const char *ty[3]={"Meas","Slack","Cnsr"};
00913 static const char *fx[2]={"Used","Fixed"};
00914
00915 if (!name) name = "";
00916 printf("TCFitV0::Print(%s) nMeas=%d(%d) nSlac=%d(%d) nCons=%d(%d)\n"
00917 ,name
00918 ,GetNMeas(),GetUMeas()
00919 ,GetNSlac(),GetUSlac()
00920 ,GetNCons(),GetUCons());
00921 for (int i=0;i<kMaxId;i++) {
00922 if (!fNams[i].Length()) continue;
00923 const double *p = 0;
00924 if (i< 5) { p = mTkFit[0].Arr()+i ;}
00925 else if (i<15) { p = mTkFit[1].Arr()+i -10;}
00926 if (!p) continue;
00927 printf("%2d - %s\t",i,GetNam(i));
00928 printf(" %s.%s ",ty[GetType(i)],fx[IsFixed(i)]);
00929
00930 printf(" %g \n",*p);
00931 }
00932 TCFitData::Print(name);
00933 }
00934
00935 #include "TRandom.h"
00936 #include "TCanvas.h"
00937 #include "TH1F.h"
00938 #include "TSystem.h"
00939
00940 static double Pdk(double ma,double mb,double mc)
00941 {
00942 double ma2 = ma*ma,ma4=ma2*ma2;
00943 double mb2 = mb*mb,mb4=mb2*mb2;
00944 double mc2 = mc*mc,mc4=mc2*mc2;
00945 double ans= (ma4+mb4+mc4-2*ma2*mb2-2*ma2*mc2-2*mb2*mc2);
00946 if (ans<0.) ans = 0;
00947 ans= sqrt(ans)/(2.*ma);
00948 return ans;
00949 }
00950
00951 void TCFitV0::Test(int mode)
00952 {
00953 static const double D0Mass=1.8646,PiMass=.13956995,KMass=.493677;
00954 static const double D0Tau=0.415e-12,C=299792458e2;
00955 static const double DX=0.003;
00956 static const double D0P=1
00957 ,D0Beta = D0P/sqrt(D0P*D0P+D0Mass*D0Mass)
00958 ,D0Gama = 1./sqrt(1.-D0Beta*D0Beta)
00959 ,D0Len = D0Tau*C*D0Gama*D0Beta;
00960 static const double HZ= (0.000299792458 * 4.98478);
00961
00962 static const int NHISTS = 4;
00963 static TCanvas* myCanvas=0;
00964 static TH1F *hh[]={0,0,0,0,0,0};
00965 static const char *hNams[]= {"dL","dL/L", "Xi2","Mass-D0"};
00966 static const double LOW[] = {-D0Len,-3, 0.,-1.5};
00967 static const double UPP[] = {-D0Len, 3,10., 1.5};
00968
00969
00970
00971 if(!myCanvas) myCanvas=new TCanvas("TCircleFitter_Test","",600,800);
00972 myCanvas->Clear();
00973 myCanvas->Divide(1,NHISTS);
00974 for (int i=0;i<NHISTS;i++) {
00975 delete hh[i]; hh[i]= new TH1F(hNams[i],hNams[i],100,LOW[i],UPP[i]);
00976 myCanvas->cd(i+1); hh[i]->Draw();
00977 }
00978
00979 TLorentzVector p4[3];
00980 TVector3 dcayPos,d0Mom;
00981 THelixTrack Khlx,Pihlx;
00982 double Q = Pdk(D0Mass,PiMass,KMass);
00983
00984 TCFitV0 dat;
00985 TCFit tc("TestD0",&dat);
00986 dat.mTkBas[0].SetHz(1.);
00987 dat.mTkBas[1].SetHz(1.);
00988
00989 for (int iev=0; iev <500; iev++) {
00990 tc.Reset(); dat.Reset();
00991 double dist = gRandom->Exp(D0Len);
00992 double wk[9];
00993 gRandom->Sphere(wk[0],wk[1],wk[2],dist);
00994 dcayPos.SetXYZ(wk[0],wk[1],wk[2]);
00995 d0Mom = dcayPos.Unit()*D0P;
00996 p4[2].SetVectM( d0Mom,D0Mass);
00997 gRandom->Sphere(wk[0],wk[1],wk[2],Q);
00998 p4[0].SetXYZM(-wk[0],-wk[1],-wk[2],KMass );
00999 p4[1].SetXYZM( wk[0], wk[1], wk[2],PiMass);
01000 TVector3 boost = p4[2].BoostVector();
01001 dat.mMas = D0Mass;
01002 for (int ip=0;ip<2;ip++) {
01003 p4[ip].Boost(boost);
01004 p4[ip].GetXYZT(wk+3);
01005 dcayPos.GetXYZ (wk+0);
01006 double ptin = 1./p4[ip].Pt();
01007 if (ip) ptin=-ptin;
01008
01009 THelixTrack th(wk+0,wk+3,ptin*HZ);
01010 th.Backward();
01011 double s =th.Path(0.,0.);
01012 th.Move(s);th.Backward();
01013 th.Eval(0.,wk,wk+3);
01014 TVector3 pos(wk),dir(wk+3);
01015 dat.mTkBas[ip].Reset();
01016 dat.mTkBas[ip].Set(pos,dir,ptin);
01017 dat.mTkBas[ip].mass= (ip)? PiMass:KMass;
01018 dat.mTkBas[ip].Update();
01019 dat.mTEBas[ip].Set(0,0,pow(DX,2));
01020 dat.mTEBas[ip].Set(1,1,pow(DX,2));
01021 dat.mTkBas[ip].Rand(dat.mTEBas[ip]);
01022 }
01023 dat.Ready();
01024 if (mode==1) dat.FixPar(kCNRJ);
01025 if (tc.Fit()) continue;
01026 double dL = dat.GetPar(kLEN_2)-dist;
01027 double eL = sqrt(dat.ErMx(kLEN_2,kLEN_2));
01028 hh[0]->Fill(dL/eL);
01029 hh[1]->Fill(dL/dist);
01030 hh[2]->Fill(dat.GetFcn()/dat.GetNDF());
01031 TLorentzVector D0V = dat.mTkFit[0].P4()+dat.mTkFit[1].P4();
01032 hh[3]->Fill(D0V.M()-D0Mass);
01033
01034 }
01035 myCanvas->Modified();
01036 myCanvas->Update();
01037
01038 }
01039
