00001
00002 #include <stdio.h>
00003 #include <iostream>
00004 #include <fstream>
00005 #include "TSystem.h"
00006 #include "TMath.h"
00007 #include "TBenchmark.h"
00008 #include "TFile.h"
00009 #include "TTree.h"
00010 #include "TH1.h"
00011 #include "TF1.h"
00012 #include "TCanvas.h"
00013 #include "TGraph.h"
00014 #include "TCernLib.h"
00015 #include "TMatrixD.h"
00016 #include "TRandom.h"
00017 #include "TVectorD.h"
00018 #include "TVector2.h"
00019 #include "TVector3.h"
00020 #include "StarRoot/TTreeIter.h"
00021 #include "TTable.h"
00022 #include "TInterpreter.h"
00023 #include "TPolinom.h"
00024 #include "THelixTrack.h"
00025 #include <vector>
00026
00027 #define SQ(X) ((X)*(X))
00028 #define PW2(X) ((X)*(X))
00029 #define PW3(X) ((X)*(X)*(X))
00030
00031 enum {kNOBAGR=0,kNDETS=4,kMINHITS=50};
00032 static const double WINDOW_NSTD=3;
00033
00034 static const double kAGREE=1e-7,kSMALL=1e-9,kBIG=1.,kDAMPR=0.1;
00035 static const int MinErr[4][2] = {{200,200},{200,200},{30,500},{20,20}};
00036 static char gTimeStamp[16]={0};
00037 class MyPull;
00038 class HitPars_t;
00039 int fiterr(const char *opt);
00040 double Fit(HitPars_t &pout);
00041 void AveRes();
00042 double AveErr();
00043 int kind(double xl);
00044 void DbInit();
00045 void DbDflt();
00046 void DbEnd();
00047 void Init(HitPars_t &hitp);
00048 void HInit();
00049 void HFit();
00050 void HEnd();
00051 void CheckErr();
00052 void FillPulls(int befAft);
00053 double newPull(int iYZ,const MyPull& myRes,const HitPars_t &pars);
00054 TBranch *tb=0;
00055 TFile *pulFile =0;
00056 TTree *pulTree=0;
00057 TH1F *hz=0;
00058 TH1F *hx=0;
00059 TH1F *hh[100];
00060 TCanvas *C[10];
00061
00062 class myTCL {
00063 public:
00064 static double vmaxa (const double *a, int na);
00065 static double vmaxa (const TVectorD &a);
00066 static void vfill (double *a,double f,int na);
00067 static void mxmlrt(const TMatrixD &A,const TMatrixD &B,TMatrixD &X);
00068 static void mxmlrtS(const TMatrixD &A,const TMatrixD &B,TMatrixD &X);
00069 static void mxmlrtS(const double *A,const double *B,double *X,int nra,int nca);
00070 static void eigen2(const double err[3], double lam[2], double eig[2][2]);
00071 static double simpson(const double *F,double A,double B,int NP);
00072 static double vasum(const double *a, int na);
00073 static double vasum(const TVectorD &a);
00074 static int SqProgSimple( TVectorD &x
00075 ,const TVectorD &g,const TMatrixD &G
00076 ,const TVectorD &Min
00077 ,const TVectorD &Max,int iAkt);
00078 };
00079
00080
00081 class HitAccr {
00082 public:
00083 int iDet;
00084 double A[2][3];
00085 double Pull[2];
00086 double PredErr[2];
00087 };
00088
00089
00090 class MyPull {
00091 public:
00092 float x_g() const { return grf*cos(gpf);}
00093 float y_g() const { return grf*sin(gpf);}
00094 float z_g() const { return gzf ;}
00095
00096 float xhg() const { return xyz[0]*cos(ang) - xyz[1]*sin(ang);}
00097 float yhg() const { return xyz[0]*sin(ang) + xyz[1]*cos(ang);}
00098 float zhg() const { return xyz[2] ;}
00099 public:
00100 int trk;
00101 float xyz[3];
00102 float psi;
00103 float dip;
00104 float pt;
00105 float ypul;
00106 float zpul;
00107 float yfit;
00108 float zfit;
00109 float pye;
00110 float pze;
00111 float uyy;
00112 float uzz;
00113 float fye;
00114 float fze;
00115 float hye;
00116 float hze;
00117 float curv;
00118 float dens;
00119 float grf;
00120 float gpf;
00121 float gzf;
00122 float ang;
00123 };
00124
00125 class HitPars_t {
00126 public:
00127 HitPars_t();
00128 HitPars_t(const HitPars_t &fr);
00129 const double &Err(int idx) const {return mDrr[idx];}
00130 double &Err(int idx) {return mDrr[idx];}
00131 double Err( int iYZ,const HitAccr &accr) const;
00132 double Err(int iDet,int iYZ,const double A[3]) const;
00133 const double &operator[](int i) const;
00134 double &operator[](int i);
00135 HitPars_t &operator= (const HitPars_t &fr);
00136 HitPars_t &operator* (double f);
00137 HitPars_t &operator+=(const double *f);
00138 HitPars_t &operator= (double f);
00139 HitPars_t &operator+=(const TVectorD &add);
00140 int NPars() const {return mNPars ;}
00141 int Len(int iDet,int iYZ=0) const {return mLen[iDet][iYZ];}
00142 int Lim(int i) const ;
00143 const double &Min(int i) const {return mMin[i];}
00144 double &Min(int i) {return mMin[i];}
00145 const double &Max(int i) const {return mMax[i];}
00146 double &Max(int i) {return mMax[i];}
00147
00148 const double &Min(int iDet,int iYZ,int i) const {return mMin[IPar(iDet,iYZ)+i];}
00149 const double &Max(int iDet,int iYZ,int i) const {return mMax[IPar(iDet,iYZ)+i];}
00150 double &Min(int iDet,int iYZ,int i) {return mMin[IPar(iDet,iYZ)+i];}
00151 double &Max(int iDet,int iYZ,int i) {return mMax[IPar(iDet,iYZ)+i];}
00152
00153 int IPar(int iDet,int iYZ,int *npars=0) const;
00154 void Set(int iDet,int iYZ,int nini,const double *ini);
00155 void Limit();
00156 double Deriv(int npt, const MyPull *pt,TVectorD &Di,TMatrixD &Dij) const;
00157 double DERIV(int npt, const MyPull *pt,TVectorD &Di,TMatrixD &Dij,int maxTrk=9999999);
00158 int Test(int npt, const MyPull *pt) const;
00159 void Print(const HitPars_t *init=0) const;
00160 double Diff (const HitPars_t &init) const;
00161 static void Prep(int npt, const MyPull *pt,TVectorD &Y,TVectorD &Z
00162 ,TVectorD &S,TVectorD &cos2Psi);
00163 static int Test();
00164 static void HitCond(const MyPull& myRes,HitAccr &acc);
00165 static void Show(int npt,const MyPull *pt);
00166 static double Dens(double rxy,int ntk);
00167 static double Err(const double Pars[3],int nPars,const double A[3]);
00168 static void myDers( double fake, double wy, double wz, double dens
00169 , double g[2][3],double G[2][3][3]);
00170 static double myFake( double fake, double wy, double wz, double dens
00171 , double Cd[3],double Cdd[3][3]);
00172 public:
00173 int mNDets;
00174 int mNPars;
00175 int mNTrks;
00176 int mLen[kNDETS][2];
00177 double *mPars[kNDETS][2];
00178 double *mErrs[kNDETS][2];
00179 double mDat[1+kNDETS*2*3];
00180 double mDrr[1+kNDETS*2*3];
00181 double mMin[1+kNDETS*2*3];
00182 double mMax[1+kNDETS*2*3];
00183
00184 private:
00185 static void myTrans( double fake,double wy, double wz, double dens, TMatrixD *mx);
00186 };
00187 HitPars_t operator+(const HitPars_t &a,const double *add);
00188 HitPars_t operator+(const HitPars_t &a,const TVectorD &add);
00189
00190
00191
00192
00193
00194 class Poli2 {
00195 public:
00196 Poli2(int npw=1);
00197 Poli2(int npw,int n,const double *X,const double *Y,const double *W);
00198 void Init();
00199 void Clear();
00200 void Add(double x, double y, double w);
00201 double Xi2() const {return fXi2;}
00202 double Xi2(int ipt) const;
00203 double EvalXi2() const;
00204 int Pw() const {return fPw ;}
00205 int NPt() const {return fN ;}
00206 double Fun( double x ) const;
00207 double dXi2dW( int ipt ) const;
00208 double d2Xi2d2W( int ipt,int jpt ) const;
00209 double Deriv( TVectorD &Di, TMatrixD &Dij ) const;
00210 void TestIt() const;
00211 void MyTest() const;
00212 static void Test();
00213 private:
00214 int fPw;
00215 char fBeg[1];
00216 int fN;
00217 double fX[100];
00218 double fY[100];
00219 double fW[100];
00220 double fXi2;
00221 double fX0,fY0;
00222 char fEnd[1];
00223 TVectorD fP;
00224 TVectorD fB;
00225 TMatrixD fA;
00226 TMatrixD fAi;
00227 };
00228
00229
00230
00231 #if !defined(__MAKECINT__)
00232
00233
00234
00235
00236
00237 std::vector<MyPull> MyVect;
00238 double aveRes[4][6],aveTrk[4][2][3];
00239 int numRes[4];
00240 double pSTI[8][3] = {{0.000421985 ,0.00124327 ,0.0257111 }
00241 ,{0.000402954 ,0.00346896 ,0.0377259 }
00242 ,{0. ,0.0009366 ,0.0004967 }
00243 ,{0.00018648 ,0.00507244 ,0.002654 }
00244 ,{0.0030*0.0030 ,0.0 ,0.0 }
00245 ,{0.0700*0.0700 ,0.0 ,0.0 }
00246 ,{0.0080*0.0080 ,0.0 ,0.0 }
00247 ,{0.0080*0.0080 ,0.0 ,0.0 }};
00248 HitPars_t HitErr;
00249
00250 static const char *DETS[]={"OutY","OutZ","InnY","InnZ"
00251 ,"SsdY","SsdZ","SvtY","SvtZ"};
00252 static const char *DETZ[]={"Out" ,"Inn","Ssd","Svt"};
00253 int FitOk[4]={0,0,0,0};
00254 static char dbFile[4][100] = {
00255 "StarDb/Calibrations/tracker/tpcOuterHitError.20050101.235959.C",
00256 "StarDb/Calibrations/tracker/tpcInnerHitError.20050101.235959.C",
00257 "StarDb/Calibrations/tracker/ssdHitError.20050101.235959.C" ,
00258 "StarDb/Calibrations/tracker/svtHitError.20050101.235959.C" };
00259
00260 static TTable *dbTab[4];
00261 void myBreak(int kase) {
00262 static int myKase=-1946;
00263 if (kase != myKase) return;
00264 printf("BOT OHO\n");
00265 }
00266
00267 int fiterr(const char *opt)
00268 {
00269 if (!opt) opt = "H";
00270 int optH = strstr(opt,"H")!=0;
00271 int optU = strstr(opt,"U")!=0;
00272 int optT = strstr(opt,"T")!=0;
00273 memcpy(gTimeStamp,strstr(opt,"20"),15);
00274
00275 DbInit();
00276 if (optH) HInit();
00277 TTreeIter th(pulTree);
00278 th.AddFile("pulls.root");
00279
00280
00281 const int &run = th("mRun");
00282 const int &evt = th("mEvt");
00283
00284 const int &mTrksG = th("mTrksG");
00285 const int &nGlobs = th("mHitsG");
00286 const float *&mCurv = th("mHitsG.mCurv");
00287 const float *&mPt = th("mHitsG.mPt");
00288 const short *&mTrackNumber = th("mHitsG.mTrackNumber");
00289 const float *&mNormalRefAngle = th("mHitsG.mNormalRefAngle");
00290 const UChar_t *&nTpcHits = th("mHitsG.nTpcHits");
00291 const UChar_t *&nSsdHits = th("mHitsG.nSsdHits");
00292 const UChar_t *&nSvtHits = th("mHitsG.nSvtHits");
00293 const float *&lYPul = th("mHitsG.lYPul");
00294 const float *&lZPul = th("mHitsG.lZPul");
00295 const float *&lXHit = th("mHitsG.lXHit");
00296 const float *&lYHit = th("mHitsG.lYHit");
00297 const float *&lZHit = th("mHitsG.lZHit");
00298 const float *&lYFit = th("mHitsG.lYFit");
00299 const float *&lZFit = th("mHitsG.lZFit");
00300 const float *&lPsi = th("mHitsG.lPsi");
00301 const float *&lDip = th("mHitsG.lDip");
00302 const float *&lYFitErr = th("mHitsG.lYFitErr");
00303 const float *&lZFitErr = th("mHitsG.lZFitErr");
00304 const float *&lYHitErr = th("mHitsG.lYHitErr");
00305 const float *&lZHitErr = th("mHitsG.lZHitErr");
00306 const float *&lYPulErr = th("mHitsG.lYPulErr");
00307 const float *&lZPulErr = th("mHitsG.lZPulErr");
00308 const float *&gRFit = th("mHitsG.gRFit");
00309 const float *&gPFit = th("mHitsG.gPFit");
00310 const float *&gZFit = th("mHitsG.gZFit");
00311
00312 printf("*lYPul = %p\n",lYPul);
00313 printf("&run=%p &evt=%p\n",&run,&evt);
00314 MyPull mp;
00315 MyVect.clear();
00316 int nTpc=0,nPre=0,iTk=-1,iTkSkip=-1;
00317 while (th.Next())
00318 {
00319
00320 float rxy=0;
00321 for (int ih=0;ih<nGlobs;ih++) {
00322 if (mTrackNumber[ih]==iTkSkip) continue;
00323 float rxy00 = rxy; rxy=lXHit[ih];
00324 if (mTrackNumber[ih]!=iTk || rxy<rxy00) {
00325 iTk = mTrackNumber[ih];
00326 iTkSkip = iTk;
00327 if (nTpcHits[ih] <25) continue;
00328 if (fabs(mCurv[ih])>1./200) continue;
00329 if (fabs(mPt[ih]/cos(lDip[ih]))<0.5) continue;
00330 int nSS = nSsdHits[ih]+nSvtHits[ih];
00331 if (nSS && nSS<2) continue;
00332 if (!nSS && nPre && nTpc>100 && nTpc>3*nPre) continue;
00333 nTpc += (nSS==0);
00334 nPre += (nSS!=0);
00335 iTkSkip = -1;
00336 }
00337
00338 memset(&mp,0,sizeof(mp));
00339 mp.trk = mTrackNumber[ih];
00340 mp.pye = lYPulErr[ih];
00341 mp.pze = lZPulErr[ih];
00342 float uyy = (lYPulErr[ih]-lYHitErr[ih])*(lYHitErr[ih]+lYPulErr[ih]);
00343
00344
00345 float uzz = (lZPulErr[ih]-lZHitErr[ih])*(lZHitErr[ih]+lZPulErr[ih]);
00346
00347
00348 mp.uyy = uyy;
00349 mp.uzz = uzz;
00350 mp.fye = lYFitErr[ih];
00351 mp.fze = lZFitErr[ih];
00352 mp.hye = lYHitErr[ih];
00353 mp.hze = lZHitErr[ih];
00354 mp.xyz[0] = lXHit[ih];
00355 if (mp.xyz[0]<4) continue;
00356 mp.xyz[1] = lYHit[ih];
00357 mp.xyz[2] = lZHit[ih];
00358 mp.psi = lPsi[ih] ;
00359 mp.dip = lDip[ih] ;
00360 mp.ypul = lYPul[ih];
00361
00362 mp.zpul = lZPul[ih];
00363
00364 mp.yfit = lYFit[ih];
00365 mp.zfit = lZFit[ih];
00366 mp.dens = HitPars_t::Dens(mp.xyz[0],mTrksG);
00367 mp.grf = gRFit[ih];
00368 mp.gpf = gPFit[ih];
00369 mp.gzf = gZFit[ih];
00370 mp.curv = mCurv[ih];
00371 mp.ang = mNormalRefAngle[ih];
00372
00373 MyVect.push_back(mp);
00374
00375 }
00376
00377 }
00378 printf("fiterr:: %d hits accepted\n\n",MyVect.size());
00379 printf("fiterr:: %d Tpc and %d Svt/Ssd tracks accepted\n\n",nTpc,nPre);
00380 AveRes();
00381 HFit();
00382 DbDflt();
00383 Init(HitErr);
00384 if (optH) CheckErr();
00385 if (optH) FillPulls(0);
00386 if (optT) return HitErr.Test(MyVect.size(),&(MyVect[0]));
00387
00388 double maxPct = Fit(HitErr);
00389 maxPct= AveErr();
00390 if (optH) FillPulls(1);
00391 if (optU) DbEnd();
00392 if (optH) HEnd();
00393
00394 return (maxPct< 3) ? 99:0;
00395 }
00396
00397 Poli2::Poli2(int npw):fP(npw+1),fB(npw+1),fA(npw+1,npw+1),fAi(npw+1,npw+1)
00398 {
00399 fPw = npw;
00400 Clear();
00401
00402 }
00403
00404 void Poli2::Clear()
00405 {
00406 memset(fBeg,0,fEnd-fBeg+1);
00407 }
00408
00409 Poli2::Poli2(int npw,int N,const double *X,const double *Y,const double *W)
00410 :fP(npw+1),fB(npw+1),fA(npw+1,npw+1),fAi(npw+1,npw+1)
00411 {
00412 fPw = npw;
00413 Clear();
00414 fN = N;
00415 int nb = N*sizeof(double);
00416 memcpy(fX,X,nb);
00417 memcpy(fY,Y,nb);
00418 memcpy(fW,W,nb);
00419 }
00420
00421 void Poli2::Add(double x, double y, double w)
00422 {
00423 fX[fN]=x; fY[fN]=y, fW[fN]=w; fN++;
00424 assert(fN<=100);
00425 }
00426
00427 void Poli2::Init()
00428 {
00429 double Wt=0,Wx=0,Wy=0;
00430 for (int i=0;i<fN;i++) { Wt += fW[i]; Wx += fW[i]*fX[i];Wy += fW[i]*fY[i];}
00431 fX0 = Wx/Wt; fY0 = Wy/Wt;
00432 for (int i=0;i<fN;i++) { fX[i]-=fX0; fY[i]-=fY0;}
00433
00434 double A[3][3]={{0}},B[3]={0};
00435 double x[3]={1};
00436 fXi2=0;
00437 for (int i=0;i<fN;i++) {
00438 double w = fW[i];
00439 x[1] = fX[i]; x[2]=x[1]*x[1];
00440 double y = fY[i],y2=y*y;
00441 for (int j=0;j<=fPw;j++) { B[j] +=w*x[j]*y;
00442 for (int k=0;k<=j ;k++) { A[j][k] +=w*x[j]*x[k];}}
00443 fXi2 +=w*y2;
00444 }
00445 for (int j=0;j<=fPw;j++){ fB[j] = B[j];
00446 for (int k=0;k<=j ;k++){ fA(j,k) = A[j][k]; fA(k,j)=A[j][k];}}
00447
00448 double det=0;
00449 fAi=fA;fAi.InvertFast(&det);
00450 fP=fAi*fB;
00451 fXi2 -= (fB*fP);
00452
00453 }
00454
00455 double Poli2::Fun( double x ) const
00456 {
00457 x -= fX0;
00458 double xx[3]={1,x,x*x};
00459 TVectorD XX(fPw+1,xx);
00460 return fP*XX + fY0;
00461 }
00462
00463 double Poli2::Xi2(int ipt) const
00464 {
00465 double dy = Fun(fX[ipt]+fX0) - (fY[ipt]+fY0);
00466 return dy*dy*fW[ipt];
00467 }
00468
00469 double Poli2::EvalXi2() const
00470 {
00471 double sum=0;
00472 for (int ipt=0;ipt<fN;ipt++) {sum+=Xi2(ipt);}
00473 return sum;
00474 }
00475
00476 double Poli2::dXi2dW( int ipt ) const
00477 {
00478
00479
00480
00481
00482 double der = fY[ipt]*fY[ipt];
00483 TVectorD dB(fPw+1);
00484 TMatrixD dA(fPw+1,fPw+1);
00485 double x[3]={1,fX[ipt],fX[ipt]*fX[ipt]};
00486 for (int j=0;j<=fPw;j++) { dB[j] = x[j]*fY[ipt];
00487 for (int k=0;k<=fPw;k++) { dA(j,k) = x[j]*x[k]; }}
00488 der -= (2*(dB*(fAi*fB))-(fP*(dA*fP)));
00489 return der;
00490 }
00491
00492
00493 double Poli2::d2Xi2d2W( int ipt,int jpt ) const
00494 {
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00505
00506 double xi[3]={1,fX[ipt],fX[ipt]*fX[ipt]};
00507 double xj[3]={1,fX[jpt],fX[jpt]*fX[jpt]};
00508
00509 TVectorD diB(fPw+1) ,djB(fPw+1);
00510 TMatrixD diA(fPw+1,fPw+1),djA(fPw+1,fPw+1);
00511 for (int j=0;j<=fPw;j++) { diB[j]=fY[ipt]*xi[j]; djB[j]=fY[jpt]*xj[j];
00512 for (int k=0;k<=fPw;k++) { diA(j,k) = xi[j]*xi[k]; djA(j,k) = xj[j]*xj[k];}}
00513 double der = -2*(-((fAi*diB)*(djA*fP))
00514 +(diB*(fAi*djB))
00515 -((fAi*djB)*(diA*fP))
00516 +(fP*(djA*(fAi*(diA*fP)))));
00517 return der;
00518 }
00519
00520 double Poli2::Deriv( TVectorD &Di, TMatrixD &Dij ) const
00521 {
00522 Di.ResizeTo(fN);
00523 Dij.ResizeTo(fN,fN);
00524 TVectorD diB(fPw+1) ,djB(fPw+1);
00525 TMatrixD diA(fPw+1,fPw+1),djA(fPw+1,fPw+1);
00526 for (int ipt=0;ipt<fN;ipt++) {
00527 double der = fY[ipt]*fY[ipt];
00528 double xi[3]={1,fX[ipt],fX[ipt]*fX[ipt]};
00529 for (int k=0;k<=fPw;k++) { diB[k] = fY[ipt]*xi[k];
00530 for (int l=0;l<=fPw;l++) { diA(k,l) = xi[k]*xi[l];}}
00531 der -= (2*(diB*(fAi*fB))-(fP*(diA*fP)));
00532 Di[ipt] = der;
00533 for (int jpt=0;jpt<=ipt;jpt++) {
00534 double xj[3]={1,fX[jpt],fX[jpt]*fX[jpt]};
00535 for (int k=0;k<=fPw;k++) { djB[k] = fY[jpt]*xj[k];
00536 for (int l=0;l<=fPw;l++) { djA(k,l) = xj[k]*xj[l];}}
00537 der = -2*(-((fAi*diB)*(djA*fP))
00538 +(diB*(fAi*djB))
00539 -((fAi*djB)*(diA*fP))
00540 +(fP*(djA*(fAi*(diA*fP)))));
00541 Dij[ipt][jpt]=der;
00542 Dij[jpt][ipt]=der;
00543 }
00544 }
00545 return Xi2();
00546 }
00547
00548
00549 void Poli2::MyTest() const
00550 {
00551 printf("Print B\n");
00552 fB.Print();
00553 printf("Print A\n");
00554 fA.Print();
00555 printf("Print Ai\n");
00556 fAi.Print();
00557 printf("Print A*Ai\n");
00558 TMatrixD tst(fPw+1,fPw+1);
00559 tst = fA*fAi;
00560 tst.Print();
00561 }
00562
00563 void Poli2::Test()
00564 {
00565 int npw=2;
00566 double X[20],Y[20],W[20];
00567 for (int i=0;i<20;i++) {
00568 X[i]=i;
00569 Y[i]= 3+X[i]*(.02+.03*X[i]);
00570 W[i]= 1-0.01*i;
00571 Y[i]+=gRandom->Gaus(0,sqrt(1./W[i]));
00572 }
00573 Poli2 pp(npw,20,X,Y,W);
00574 pp.Init();
00575
00576 double Xi2 = pp.Xi2();
00577 double Xi2Eva = pp.Xi2();
00578 printf ("Xi2 = %g == %g\n",Xi2,Xi2Eva);
00579 for (int i=0; i<20; i++) {
00580 printf(" %g %g : %g\n",X[i],Y[i],pp.Fun(X[i]));
00581 }
00582
00583 printf("\n check d/dWi\n");
00584 double delta = 1e-3;
00585 for (int k=0;k<20;k++) {
00586 Poli2 ppp(npw);
00587 for (int i=0;i<20;i++) {
00588 double w = W[i];
00589 if (i==k) w+=delta;
00590 ppp.Add(X[i],Y[i],w);
00591 }
00592 ppp.Init();
00593 double ana = pp.dXi2dW(k);
00594 double num = (ppp.Xi2()-pp.Xi2())/delta;
00595 double dif = 2*(num-ana)/(fabs(num+ana)+3e-33);
00596 if (fabs(dif)<0.01) continue;
00597 printf ("dXi2dW(%2d) \tana=%g \tnum = %g \tdif=%g\n",k,ana,num,dif);
00598 }
00599
00600 printf("\n check d2/dWi/DWj\n");
00601 for (int i=0;i<20;i++) {
00602 double sav = W[i]; W[i]=sav+delta;
00603 Poli2 ppp(npw,20,X,Y,W);
00604 W[i]=sav;
00605 ppp.Init();
00606 for (int j=0;j<20;j++) {
00607 double ana = pp.d2Xi2d2W(i,j);
00608 double num = (ppp.dXi2dW(j)-pp.dXi2dW(j))/delta;
00609 double dif = 2*(num-ana)/(fabs(num+ana)+3e-33);
00610 if (fabs(dif)<0.01) continue;
00611 printf ("d2Xi2d2W(%2d,%2d) \tana=%g \tnum = %g \tdif=%g\n",i,j,ana,num,dif);
00612 } }
00613 TVectorD g;
00614 TMatrixD gg;
00615 pp.Deriv(g,gg);
00616 for (int i=0;i<20;i++) {
00617 double tst = pp.dXi2dW(i)-g[i];
00618 if (fabs(tst)>1e-10) printf("g[%d] %g **************\n",i,tst);
00619 for (int j=0;j<20;j++) {
00620 double tst = pp.d2Xi2d2W(i,j)-gg[i][j];
00621 if (fabs(tst)>1e-10) printf("gg[%d][%d] %g **************\n",i,j,tst);
00622 } }
00623
00624
00625 }
00626
00627 void Poli2::TestIt() const
00628 {
00629 TVectorD g,gT;
00630 TMatrixD G,GT;
00631 Poli2 pp(*this);
00632 pp.Init();
00633 double Xi2 = pp.Deriv(g,G);
00634
00635 double delta = 1e-3;
00636 for (int k=0;k<fN;k++) {
00637 double sav = pp.fW[k];
00638 pp.fW[k]=sav+delta;
00639 pp.Init();
00640 double Xi2T = pp.Deriv(gT,GT);
00641 pp.fW[k]=sav;
00642
00643 double ana = 0.5*(g[k]+gT[k]);
00644 double num = (Xi2T-Xi2)/delta;
00645 double dif = 2*(num-ana)/(fabs(num+ana)+1e-10);
00646 if (fabs(dif)>0.01)
00647 printf ("\ndXi2dW(%2d) \tana=%g \tnum = %g \tdif=%g\n\n",k,ana,num,dif);
00648
00649 for (int i=0;i<fN;i++) {
00650 ana = 0.5*(GT[k][i]+G[k][i]);
00651 num = (gT[i]-g[i])/delta;
00652 dif = 2*(num-ana)/(fabs(num+ana)+1e-10);
00653 if (fabs(dif)>0.01)
00654 printf ("d2Xi2dW2(%2d,%2d) \tana=%g \tnum = %g \tdif=%g\n",k,i,ana,num,dif);
00655 } }
00656
00657
00658 }
00659
00660 void poli2()
00661 {
00662 Poli2::Test();
00663 }
00664
00665 void Init(HitPars_t &hiterr)
00666 {
00667 hiterr[0] = 0.0;
00668 for (int iDet=0;iDet<kNDETS;iDet++) {
00669 if (numRes[iDet]<kMINHITS) continue;
00670 int n = (iDet<=1)? 3:1;
00671 hiterr.Set(iDet,0,n,pSTI[iDet*2+0]);
00672 hiterr.Set(iDet,1,n,pSTI[iDet*2+1]);
00673 hiterr.Min(iDet,0,0) = pow(MinErr[iDet][0]*1.e-4,2);
00674 hiterr.Min(iDet,1,0) = pow(MinErr[iDet][1]*1.e-4,2);
00675 }
00676 }
00677
00678
00679 class FitState_t : public HitPars_t {
00680 public:
00681 FitState_t();
00682 FitState_t(HitPars_t &hp);
00683 const HitPars_t &Pars() const { return (const HitPars_t &)(*this);}
00684 HitPars_t &Pars() { return ( HitPars_t &)(*this);}
00685 int Saddle() const {return neg<0;}
00686 int LimX(int i) const;
00687 double Der() const {return der;}
00688 double Fcn() const {return fcn;}
00689 void Deriv(const std::vector<MyPull> &MyVect);
00690 int MaxStp(TVectorD &add,int mode) const;
00691 void MakeErrs();
00692 int operator<(const FitState_t &other) const;
00693 FitState_t &operator=(const FitState_t &other);
00694 static int FixWeak( TVectorD &g, TMatrixD &G);
00695 public:
00696 TVectorD g;
00697 TMatrixD G;
00698 char myBeg[1];
00699 int npt;
00700 int ider;
00701 double fak;
00702 double fcn;
00703 double der;
00704 double neg;
00705 char myEnd[1];
00706 };
00707
00708 FitState_t::FitState_t()
00709 {
00710 memset(myBeg,0,myEnd-myBeg+1);
00711 der = 1e99;fcn = 1e99;neg =-1e99;
00712 }
00713
00714 FitState_t::FitState_t(HitPars_t &hp) :HitPars_t(hp)
00715 {
00716 memset(myBeg,0,myEnd-myBeg+1);
00717 der = 1e99;fcn = 1e99;neg =-1e99;
00718 }
00719
00720 void FitState_t::Deriv(const std::vector<MyPull> &MyVect)
00721 {
00722 npt = MyVect.size();
00723 fcn = DERIV(npt,&(MyVect[0]),g,G);
00724 fcn/=npt;
00725 der = -1; ider =-1; neg = 0;
00726 for (int i=0;i<mNPars;i++) {
00727 if (G[i][i]<neg) neg =G[i][i];
00728 if(LimX(i)) continue;
00729 if (fabs(g[i])>der) { ider = i; der = fabs(g[i]);}
00730 }
00731 der/=npt;
00732 double myMax=myTCL::vmaxa(G.GetMatrixArray(),G.GetNoElements());
00733 fak = pow(2.,-int(log(myMax)/log(2.)));
00734 G*=fak; g*=fak;
00735 if (neg>=0) {
00736 for (int i = 0; i<mNPars && neg>0; i++){
00737 for (int j = 0; j<i && neg>0; j++) {
00738 if (pow(G[i][j],2) >= G[i][i]*G[j][j]) neg = -1;
00739 } }
00740 }
00741 }
00742
00743
00744
00745 int FitState_t::operator<(const FitState_t &old) const
00746 {
00747
00748 if (Der()<0 ) return 0;
00749 if (Fcn()<old.Fcn()) return 1;
00750 if (Der()<old.Der()
00751 && Fcn()-old.Fcn()-0.1*fabs(old.Fcn())<0) return 1;
00752 return 0;
00753 }
00754
00755 FitState_t &FitState_t::operator=(const FitState_t &other)
00756 {
00757 memcpy(myBeg,other.myBeg,myEnd-myBeg+1);
00758 g.ResizeTo(other.g);
00759 G.ResizeTo(other.G);
00760 g = other.g;
00761 G = other.G;
00762 HitPars_t::operator=(other);
00763 return *this;
00764 }
00765
00766 int FitState_t::MaxStp(TVectorD &add,int mode) const
00767 {
00768 int lim=0;
00769 double fak=1;
00770 for (int i=0;i<mNPars;i++) {
00771 double maxStp = (0.01+mDat[i])*0.3;
00772 if (maxStp>fak*fabs(add[i])) continue;
00773 lim = lim*100+i+1;
00774 if (!mode) { add[i] = (add[i]<0)? -maxStp:maxStp;}
00775 else { fak = maxStp/fabs(add[i]);}
00776 }
00777 if (fak<1) add*=fak;
00778 return lim;
00779 }
00780
00781 void FitState_t::MakeErrs()
00782 {
00783
00784 TMatrixD E(G);
00785 for (int i=0;i<mNPars;i++) {
00786 if (!Lim(i)) continue;
00787 for (int j=0;j<mNPars;j++) {E(i,j) = 0; E(j,i) = 0;}; E(i,i) = 1;
00788 }
00789 double det=12345; E.InvertFast(&det);
00790 E*=fak;
00791 for (int i=0;i<mNPars;i++) {Err(i) = sqrt(E(i,i));}
00792 }
00793
00794 int FitState_t::LimX(int i) const
00795 {
00796 int lim = Lim(i);
00797 if (!lim) return 0;
00798 if (lim<0 && g[i]<0) return 0;
00799 if (lim>0 && g[i]>0) return 0;
00800 return lim;
00801 }
00802
00803 int FitState_t::FixWeak( TVectorD &g, TMatrixD &G)
00804 {
00805 int n = g.GetNrows();
00806 int nfix = 0;
00807 double ave = myTCL::vasum(g)/n;
00808 for (int i=0;i<n;i++) {
00809 if (fabs(g[i])>=ave) continue;
00810 nfix++; g[i]=0;
00811 for (int j=0;j<n;j++) {G[i][j]=0; G[j][i]=0;}
00812 G[i][i]=1;
00813 }
00814 return nfix;
00815 }
00816
00817 double Fit(HitPars_t &pout)
00818 {
00819 enum {kMAXITER=10000};
00820 static int const MAXCUT[2]={4,10};
00821 static int kount=0;
00822
00823 int ifail = 99,nPars,iAkt=0,iter,nCut,lim=0,con=0;
00824 nPars = pout.NPars();
00825 HitPars_t init(pout);
00826 init.Limit();
00827
00828 TVectorD add(nPars),g(nPars),ge(nPars),Gg(nPars);
00829 FitState_t Best(init),Now(init);
00830 double dif,difFcn=0,difDer;
00831
00832 static int idebug = 1;
00833
00834 iAkt=0;nCut=0;
00835
00836 for (iter=0;iter<kMAXITER;iter++) {
00837 kount++;
00838 Now.Deriv(MyVect);
00839 difFcn = Now.fcn-Best.fcn;
00840 difDer = Now.der-Best.der;
00841 if (idebug) {
00842 printf("Fit(%3d) \tFcn = %g(%g) \tDer(%2d)=%g(%g) \tLim=%d \tCon=%d\tAkt=%d Cut=%d Sad=%d\n"
00843 ,iter,Now.fcn,difFcn,Now.ider,Now.der,difDer,lim,con,iAkt,nCut,Now.Saddle());
00844 }
00845 if (Now.Der()<0) {ifail=1; break;
00846 } else if (Now.Der() < kAGREE) {
00847 Best=Now; ifail=0; break;
00848
00849 } else if (Now < Best ) {
00850 nCut = 0; Best = Now; iAkt = 0;
00851
00852 } else if (nCut>MAXCUT[iAkt]) {
00853 if (iAkt==0) Now=Best;
00854 Best=Now; nCut=0; iAkt=1-iAkt;
00855
00856 } else {
00857 nCut++;
00858 add*= 0.5;
00859 Now.Pars() = Best.Pars()+add;
00860 continue;
00861 }
00862
00863
00864
00865
00866 static int rnd=0; rnd++;
00867 if (!(rnd&3)) iAkt=1;
00868
00869
00870 TVectorD P0(nPars,&Best[0]);
00871 for (int jAkt=iAkt;jAkt<2;jAkt++) {
00872 iAkt = jAkt;
00873 TVectorD P1(P0);
00874 con = myTCL::SqProgSimple(P1,Now.g,Now.G
00875 ,TVectorD(nPars,&Now.Min(0))
00876 ,TVectorD(nPars,&Now.Max(0)),jAkt);
00877 if (con<0) continue;
00878 add = P1-P0;
00879 double along = -(add*Now.g);
00880 if (along <0) continue;
00881 lim=Best.MaxStp(add,1);
00882 Now.Pars() = Best.Pars()+add;
00883 break;
00884 }
00885 }
00886
00887 if (ifail==0 || ifail==99) Best.MakeErrs();
00888 pout=Best.Pars();
00889
00890 dif = pout.Diff(init);
00891 printf("\nFit: Iter=%d Fcn=%g Der=%g Dif=%6.3g%% Fail=%d\n"
00892 ,iter,Best.fcn,Best.der,dif,ifail);
00893 pout.Print(&init);
00894 return (ifail)? 1e10:dif;
00895
00896 }
00897
00898
00899 void AveRes()
00900 {
00901 memset(aveRes[0] ,0,sizeof(aveRes));
00902 memset(aveTrk[0][0],0,sizeof(aveTrk));
00903 memset(numRes ,0,sizeof(numRes));
00904 for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
00905 MyPull &myRes = MyVect[jhit];
00906 int jdx = kind(myRes.xyz[0]);
00907 aveRes[jdx][0]+= myRes.ypul*myRes.ypul;
00908 aveRes[jdx][1]+= myRes.zpul*myRes.zpul;
00909 aveRes[jdx][2]+= myRes.uyy;
00910 aveRes[jdx][3]+= myRes.uzz;
00911 aveRes[jdx][4]+= pow(myRes.xyz[1]-myRes.yfit,2);
00912 aveRes[jdx][5]+= pow(myRes.xyz[2]-myRes.zfit,2);
00913 if (hh[jdx*2+0+30]) hh[jdx*2+0+30]->Fill(myRes.ypul);
00914 if (hh[jdx*2+1+30]) hh[jdx*2+1+30]->Fill(myRes.zpul);
00915 numRes[jdx]++;
00916 HitAccr ha;
00917 HitPars_t::HitCond(myRes,ha);
00918 TCL::vadd(aveTrk[jdx][0],ha.A[0],aveTrk[jdx][0],3);
00919 TCL::vadd(aveTrk[jdx][1],ha.A[1],aveTrk[jdx][1],3);
00920
00921 }
00922
00923 int ihit=0;
00924 for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
00925 MyPull &myRes = MyVect[jhit];
00926 int jdx = kind(myRes.xyz[0]);
00927 if (numRes[jdx]<kMINHITS) continue;
00928 aveRes[jdx][0]+= myRes.ypul*myRes.ypul;
00929 if (ihit!=jhit) MyVect[ihit]=MyVect[jhit];
00930 ihit++;
00931 }
00932 MyVect.resize(ihit);
00933
00934 for(int jdx=0;jdx<4;jdx++) {
00935 if (numRes[jdx]<kMINHITS) continue;
00936 double f = 1./numRes[jdx];
00937 TCL::vscale(aveRes[jdx],f,aveRes[jdx],6);
00938 TCL::vscale(aveTrk[jdx][0],f,aveTrk[jdx][0],6);
00939 double hitYErr,hitZErr;
00940 TCL::vsub(aveRes[jdx]+0,aveRes[jdx]+2,aveRes[jdx]+2,2);
00941 hitYErr = aveRes[jdx][2];
00942 hitYErr = (hitYErr<0)? -sqrt(-hitYErr):sqrt(hitYErr);
00943 hitZErr = aveRes[jdx][3];
00944 hitZErr = (hitZErr<0)? -sqrt(-hitZErr):sqrt(hitZErr);
00945 printf("AveRes:: N=%5d \taveRes[%s][yz]=%g %g \tavePul[yz]=%g %g \thitErr(yz)=%g %g\n\n"
00946 ,numRes[jdx],DETZ[jdx]
00947 ,sqrt(aveRes[jdx][4]),sqrt(aveRes[jdx][5])
00948 ,sqrt(aveRes[jdx][0]),sqrt(aveRes[jdx][1])
00949 ,hitYErr,hitZErr);
00950
00951
00952
00953 }
00954 }
00955
00956 int kind(double x)
00957 {
00958 static const double radds[5]={120,25,16,0,0};
00959 for(int jdx=0;1;jdx++) {if (x>radds[jdx]) return jdx;}
00960 }
00961
00962 void HInit()
00963 {
00964 memset(hh,0,sizeof(hh));
00965 memset(C,0,sizeof(C));
00966
00967 static const char *NamPuls[]={
00968 "YOut.bef","YOut.aft","ZOut.bef","ZOut.aft",
00969 "YInn.bef","YInn.aft","ZInn.bef","ZInn.aft",
00970 "YSsdBef","YSsdAft","ZSsdBef","ZSsdAft",
00971 "YSvtBef","YSvtAft","ZSvtBef","ZSvtAft"};
00972 static const char *NamRes[]={
00973 "YOut.res","ZOut.res",
00974 "YInn.res","ZInn.res",
00975 "YSsd.res","ZSsd.res",
00976 "YSvt.res","ZSvt.res"};
00977
00978 for (int ih=0;ih<16;ih++) {
00979 hh[ih] = new TH1F(NamPuls[ih],NamPuls[ih],100,-3,3);}
00980 C[0] = new TCanvas("C0","",600,800);
00981 C[0]->Divide(1,8);
00982 C[1] = new TCanvas("C1","",600,800);
00983 C[1]->Divide(1,8);
00984 for (int ih=0;ih<8;ih++) {C[0]->cd(ih+1);hh[ih+0]->Draw();}
00985 for (int ih=0;ih<8;ih++) {C[1]->cd(ih+1);hh[ih+8]->Draw();}
00986
00987 hh[20]= new TH1F("ChekErr","ChekErr",100,-100.,100.);
00988 hh[21]= new TH1F("ChekPul","ChekPul",100,-100.,100.);
00989 C[2] = new TCanvas("C2","",600,800);
00990 C[2]->Divide(1,2);
00991 for (int ih=0;ih<2;ih++) {C[2]->cd(ih+1);hh[ih+20]->Draw();}
00992
00993
00994 for (int ih=0;ih<8;ih++) {
00995 hh[ih+30] = new TH1F(NamRes[ih],NamRes[ih],100,-0.3,0.3);}
00996 C[3] = new TCanvas("C3","",600,800);
00997 C[3]->Divide(1,8);
00998 for (int ih=0;ih<8;ih++) {C[3]->cd(ih+1);hh[ih+30]->Draw();}
00999
01000 }
01001
01002
01003 void FillPulls(int befAft)
01004 {
01005 MyPull myRes;
01006 for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
01007 myRes = MyVect[jhit];
01008 int jdx = kind(myRes.xyz[0]);
01009 if (befAft==0) {
01010 hh[jdx*4+0]->Fill(myRes.ypul/(myRes.pye));
01011 hh[jdx*4+2]->Fill(myRes.zpul/(myRes.pze));
01012 } else {
01013 hh[jdx*4+1]->Fill(newPull(0,myRes,HitErr));
01014 hh[jdx*4+3]->Fill(newPull(1,myRes,HitErr));
01015 }
01016 }
01017 }
01018
01019 double newPull(int iYZ,const MyPull& myRes,const HitPars_t &pars)
01020 {
01021 HitAccr accr;
01022 HitPars_t::HitCond(myRes,accr);
01023 double S = pars.Err(iYZ,accr)+accr.PredErr[iYZ];
01024 if (S<1e-6) S=1e-6;
01025 return accr.Pull[iYZ]/sqrt(S);
01026 }
01027
01028
01029 void DbInit()
01030 {
01031 gSystem->Load("libStDb_Tables.so");
01032 TString command;
01033 TTable *newdat = 0;
01034 for (int idb=0;idb<4;idb++) {
01035 memcpy(strstr(dbFile[idb],"20"),gTimeStamp,15);
01036 int ready=0;
01037 if (!gSystem->AccessPathName(dbFile[idb])) {
01038 command = ".L "; command += dbFile[idb];
01039 gInterpreter->ProcessLine(command);
01040 newdat = (TTable *) gInterpreter->Calc("CreateTable()");
01041 command.ReplaceAll(".L ",".U ");
01042 gInterpreter->ProcessLine(command);
01043 ready = 2006;
01044 } else {
01045 newdat = (TTable *)gInterpreter->Calc("new St_HitError(\"someHitError\",1)");
01046 newdat->SetUsedRows(1);
01047 }
01048 assert(newdat);
01049 dbTab[idb] = newdat;
01050 }
01051 }
01052
01053 void DbDflt()
01054 {
01055
01056
01057 for (int idb=0;idb<4;idb++)
01058 {
01059 if (numRes[idb]<kMINHITS) continue;
01060 double *d = (double *)dbTab[idb]->GetArray();
01061 if (d[0]<=0) {d[0]=aveRes[idb][0];d[3]=aveRes[idb][1];}
01062 TCL::ucopy(d+0,pSTI[idb*2+0],3);
01063 TCL::ucopy(d+3,pSTI[idb*2+1],3);
01064 }
01065 }
01066
01067 void DbEnd()
01068 {
01069 double par[2][3];
01070 for (int idb=0;idb<kNDETS;idb++)
01071 {
01072 memset(par,0,sizeof(par));
01073 if (!HitErr.mPars[idb][0]) continue;
01074 if (!HitErr.Len(idb)) continue;
01075 TString ts(dbFile[idb]);
01076 ts +=".BAK";
01077 gSystem->Rename(dbFile[idb],ts);
01078 int ny = HitErr.Len(idb,0);
01079 TCL::ucopy(HitErr.mPars[idb][0],par[0],ny);
01080 int nz = HitErr.Len(idb,1);
01081 TCL::ucopy(HitErr.mPars[idb][1],par[1],nz);
01082
01083 TCL::ucopy(par[0],(double*)dbTab[idb]->GetArray()+0,3+3);
01084
01085 std::ofstream ofs(dbFile[idb]);
01086 dbTab[idb]->SavePrimitive(ofs);
01087 }
01088 }
01089
01090
01091 void HEnd()
01092 {
01093 for (int ic=0;ic<(int)(sizeof(C)/sizeof(void*));ic++) {
01094 TCanvas *cc = C[ic]; if (!cc) continue;
01095 cc->Modified();cc->Update();
01096 }
01097 while(!gSystem->ProcessEvents()){};
01098 }
01099
01100
01101 void CheckErr()
01102 {
01103 for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
01104 MyPull &myRes = MyVect[jhit];
01105 int jdx = kind(myRes.xyz[0]); if(jdx){};
01106 if (jdx != 2) continue;
01107 float pye = myRes.pye;
01108 float hye = myRes.hye;
01109
01110 float ypul = myRes.ypul;
01111 if (fabs(ypul/pye)<2) {
01112 double tmp = (ypul*ypul-pye*pye)/(hye*hye);
01113 hh[20]->Fill(tmp);
01114 }
01115 hh[21]->Fill(pye/hye);
01116 }
01117 }
01118
01119 double AveErr()
01120 {
01121
01122
01123
01124
01125
01126
01127
01128
01129 double pctMax=0;
01130 for (int iDet=0;iDet<kNDETS;iDet++) {
01131 for (int iYZ=0;iYZ<2;iYZ++) {
01132 int idx = iYZ + 2*iDet;
01133 int nPars = HitErr.Len(iDet,iYZ);
01134 double sum=0,sum2=0,sig[3]={0};
01135 int nTot=0;
01136 for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
01137 MyPull &myRes = MyVect[jhit];
01138 int jdx = kind(myRes.xyz[0]);
01139 if (jdx != iDet) continue;
01140 nTot++;
01141 HitAccr accr;
01142 HitPars_t::HitCond(myRes,accr);
01143 double S = HitErr.Err(iYZ,accr);
01144 double delta = pow((&myRes.ypul)[iYZ],2);
01145 double P0 = delta/((&myRes.uyy)[iYZ]+pow((&myRes.hye)[iYZ],2));
01146 double P1 = delta/((&myRes.uyy)[iYZ]+HitPars_t::Err(pSTI[idx],nPars,accr.A[iYZ]));
01147 double P2 = delta/((&myRes.uyy)[iYZ]+S);
01148 sum +=S; sum2 += S*S; sig[0]+=P0; sig[1]+=P1;sig[2]+=P2;
01149 }
01150 if (!nTot) continue;
01151 sum/=nTot; sum2/=nTot; sum2-=sum*sum; if (sum2<0) sum2=0; sum2=sqrt(sum2);
01152 sum = sqrt(sum); sum2 /= 2*sum;
01153 for (int i=0;i<3;i++) {sig[i]=sqrt(sig[i]/nTot);}
01154
01155 int isig = int(10000*sum +0.5);
01156 int esig = int(10000*sum2+0.5);
01157 double newErr = sum;
01158 double oldErr = sqrt(HitPars_t::Err(pSTI[idx],nPars,aveTrk[iDet][iYZ]));
01159 double pct = 200*fabs(newErr-oldErr)/(newErr+oldErr);
01160 if (pct>pctMax) pctMax=pct;
01161 int iold = int(10000*oldErr +0.5);
01162
01163
01164 printf("AveErr(%s): <Err.old> = %4d <Err.new> %4d(+-%4dmu) Dif=%4.1g%%\t// evts=%d\n"
01165 ,DETS[idx],iold,isig,esig,pct,nTot);
01166 } }
01167 printf("AveErr(): maxPct = %4.1g%%\n\n", pctMax );
01168
01169
01170 return pctMax;
01171 }
01172
01173 void HFit()
01174 {
01175 if (!hh[30+4]) return;
01176 }
01177
01178
01179
01180
01181 HitPars_t::HitPars_t ()
01182 {
01183 memset(this,0,sizeof(HitPars_t));
01184 mNPars=1;
01185 int n = sizeof(mMin)/sizeof(*mMin);
01186 myTCL::vfill(mMin,kSMALL,n);
01187 myTCL::vfill(mMax,kBIG ,n);
01188 mMax[0] = 3;
01189 mDat[0] = 0.1;
01190 }
01191
01192 HitPars_t::HitPars_t (const HitPars_t &fr)
01193 {
01194 *this = fr;
01195 }
01196
01197 HitPars_t &HitPars_t::operator=(const HitPars_t &fr)
01198 {
01199 memcpy(this,&fr,sizeof(HitPars_t));
01200 int inc = (char*)mDat-(char*)fr.mDat;
01201 int n = sizeof(mPars)/sizeof(mPars[0][0]);
01202 char **p = (char**)mPars[0];
01203 for (int i=0;i<n;i++) {if (p[i]) p[i] += inc;}
01204 assert(mDat <=mPars[0][0]);
01205 return *this;
01206 }
01207
01208 const double &HitPars_t::operator[](int i) const
01209 { return mDat[i];}
01210
01211
01212 double &HitPars_t::operator[](int i)
01213 { return mDat[i];}
01214
01215 HitPars_t &HitPars_t::operator*(double f)
01216 { for (int i=0;i<mNPars;i++) { (*this)[i]*=f;} return *this; }
01217
01218 HitPars_t &HitPars_t::operator+=(const double *f)
01219 { for (int i=0;i<mNPars;i++) {(*this)[i]+=f[i];} return *this;}
01220
01221 HitPars_t &HitPars_t::operator=(double f)
01222 { for (int i=0;i<mNPars;i++) {(*this)[i]=f;} return *this;}
01223
01224 HitPars_t operator+(const HitPars_t &a,const double *add)
01225 {
01226 HitPars_t tmp(a); tmp+=add; return tmp;
01227 }
01228
01229 HitPars_t operator+(const HitPars_t &a,const TVectorD &add)
01230 {
01231 HitPars_t tmp(a); tmp+=add; return tmp;
01232 }
01233
01234 int HitPars_t::IPar(int iDet,int iYZ, int *npars) const
01235 {
01236 if (npars) *npars=Len(iDet,iYZ);
01237 int ans = mPars[iDet][iYZ]-mDat;
01238 assert(ans>=0 && ans < 100);
01239 return ans;
01240 }
01241
01242 int HitPars_t::Lim(int i) const
01243 {
01244 int lim = 0;
01245 if (mDat[i] <= mMin[i]*1.1) lim = -1;
01246 if (mDat[i] >= mMax[i]*0.9) lim = 1;
01247 return lim;
01248 }
01249
01250 double HitPars_t::Err(int iDet,int iYZ,const double A[3]) const
01251 {
01252 return Err(mPars[iDet][iYZ],Len(iDet,iYZ),A);
01253 }
01254
01255 double HitPars_t::Err(int iYZ,const HitAccr &accr) const
01256 {
01257 return Err(accr.iDet,iYZ,accr.A[iYZ]);
01258 }
01259
01260 HitPars_t &HitPars_t::operator+=(const TVectorD &add)
01261 {
01262 for (int i=0;i<mNPars;i++) {
01263 mDat[i]+=add[i];
01264 if (mDat[i]<mMin[i]) mDat[i]=mMin[i];
01265 if (mDat[i]>mMax[i]) mDat[i]=mMax[i];
01266 }
01267 return *this;
01268 }
01269
01270 void HitPars_t::Set(int iDet,int iYZ,int nini,const double *ini)
01271 {
01272 mLen[iDet][iYZ]=nini;
01273 if (iDet+1>mNDets) mNDets=iDet+1;;
01274 mPars[iDet][iYZ]=mDat+mNPars;
01275 mErrs[iDet][iYZ]=mDrr+mNPars;
01276 mNPars+=nini;
01277 for (int i=0;i<nini;i++) {
01278 int ipar = IPar(iDet,iYZ);
01279 assert(ipar>=0);
01280 mDat[ipar+i]=ini[i];
01281 if (mDat[ipar+i]< mMin[ipar+i]) mDat[ipar+i]= mMin[ipar+i];
01282 if (mDat[ipar+i]> mMax[ipar+i]) mDat[ipar+i]= mMax[ipar+i];
01283 }
01284 }
01285
01286
01287 double HitPars_t::Dens(double rxy,int ntk)
01288 {
01289
01290
01291
01292 return ntk/(4*3.14*rxy*rxy);
01293 }
01294
01295 void HitPars_t::HitCond(const MyPull& myRes,HitAccr &acc)
01296 {
01297 memset(acc.A[0],0,sizeof(acc.A));
01298 acc.iDet = kind(myRes.xyz[0]);
01299 for (int iyz=0;iyz<2;iyz++) {
01300 acc.A[iyz][0]=1;
01301 acc.A[iyz][1]= 0.01*(200-fabs(myRes.xyz[2]));
01302 double ang;
01303 if (!iyz) {
01304 ang = myRes.psi; acc.Pull[0]=myRes.ypul;
01305 acc.PredErr[0] = myRes.uyy;
01306 } else {
01307 ang = myRes.dip; acc.Pull[1]=myRes.zpul;
01308 acc.PredErr[1] = myRes.uzz;
01309 }
01310 double ca=cos(ang);
01311 if (ca<0.01) ca=0.01;
01312 double ca2 = ca*ca;
01313 double ta2=(1-ca2)/ca2;
01314 acc.A[iyz][1]/= ca2;
01315 acc.A[iyz][2] = ta2;
01316 }
01317 }
01318
01319 double HitPars_t::Deriv(int npt, const MyPull *pt,TVectorD &Di,TMatrixD &Dij) const
01320 {
01321
01322
01323
01324
01325
01326
01327
01328
01329
01330
01331 Di.ResizeTo(mNPars);
01332 Dij.ResizeTo(mNPars,mNPars);
01333 int npt21 = npt*2+1;
01334
01335 HitAccr ha;
01336 Poli2 ply(2);
01337 Poli2 plz(1);
01338 TVectorD wy(npt),wz(npt),cos2Psi(npt21);
01339 TMatrixD toPars(mNPars,npt21);
01340 toPars[0][0]=1;
01341 cos2Psi=1.;
01342 TPoliFitter pfy(2),pfz(1);
01343 TCircleFitter cf;
01344 #define NEWPOLIFIT
01345 #ifndef NEWPOLIFIT
01346 double len = 0,oldCurv,oldXyz[3],nowXyz[3];
01347 assert(pt[0].xyz[0]<pt[npt-1].xyz[0]);
01348 for (int ipt=0;ipt<npt;ipt++) {
01349 HitCond(pt[ipt],ha);
01350 int iDet = ha.iDet;
01351 if (!mPars[iDet][0]) continue;
01352 double nowCurv=pt[ipt].curv;
01353 nowXyz[0] = pt[ipt].grf*cos(pt[ipt].gpf);
01354 nowXyz[1] = pt[ipt].grf*sin(pt[ipt].gpf);
01355 if (ipt) {
01356 double dl = sqrt(SQ(nowXyz[0]-oldXyz[0])+SQ(nowXyz[1]-oldXyz[1]));
01357 double curv=0.5*fabs(nowCurv+oldCurv);
01358 if (dl*curv<1e-3) dl*=(1.+SQ(dl*curv)/24);
01359 else dl = 2*asin(0.5*dl*curv)/curv;
01360 len+=dl;
01361 }
01362
01363 double dy = pt[ipt].xyz[1]-pt[ipt].yfit;
01364 double cosPsi = cos(pt[ipt].psi);
01365 cos2Psi[ipt+1] = cosPsi*cosPsi;
01366 dy *= cosPsi;
01367 double dz = pt[ipt].xyz[2]-pt[ipt].zfit;
01368 wy[ipt] = (1./Err(0,ha))/cos2Psi[ipt+1];
01369 wz[ipt] = (1./Err(1,ha));
01370 double fake = myFake( mDat[0],wy[ipt],wz[ipt],pt[ipt].dens,0,0);
01371 ply.Add(len,dy,wy[ipt]*fake);
01372 plz.Add(len,dz,wz[ipt]*fake);
01373
01374
01375
01376
01377 int n,ipar;
01378 for (int jk=0;jk<2;jk++) {
01379 int ip = 1 +ipt + jk*npt;
01380 ipar = IPar(iDet,jk,&n);
01381 assert(ipar>=0);
01382 for (int in=0;in<n;in++) {toPars[ipar+in][ip]+=ha.A[jk][in];}
01383 }
01384 memcpy(oldXyz,nowXyz,sizeof(oldXyz));
01385 oldCurv=nowCurv;
01386 }
01387 #endif //0
01388 #ifdef NEWPOLIFIT
01389 TVectorD Y,Z,S;
01390 Prep(npt,pt,Y,Z,S,cos2Psi);
01391 for (int ipt=0;ipt<npt;ipt++) {
01392 HitCond(pt[ipt],ha);
01393 wy[ipt] = (1./Err(0,ha))/cos2Psi[ipt+1];
01394 wz[ipt] = (1./Err(1,ha));
01395 double fake = myFake( mDat[0],wy[ipt],wz[ipt],pt[ipt].dens,0,0);
01396 ply.Add(S[ipt],Y[ipt],wy[ipt]*fake);
01397 plz.Add(S[ipt],Z[ipt],wz[ipt]*fake);
01398 pfy.Add(S[ipt],Y[ipt],1./(wy[ipt]*fake));
01399 pfz.Add(S[ipt],Z[ipt],1./(wz[ipt]*fake));
01400 double emx[3]={1./(wy[ipt]*fake),0,1./(wy[ipt]*fake)};
01401 cf.Add(S[ipt],Y[ipt],emx);
01402 int n,ipar,iDet = ha.iDet;
01403 static int myCall=0; myCall++;
01404 for (int jk=0;jk<2;jk++) {
01405 int ip = 1 +ipt + jk*npt;
01406 ipar = IPar(iDet,jk,&n);
01407 assert(ipar>=0);
01408 for (int in=0;in<n;in++) {toPars[ipar+in][ip]+=ha.A[jk][in];}
01409 }
01410 }
01411 #endif //0
01412 ply.Init();
01413 plz.Init();
01414
01415 static int testIt=1;
01416 if (testIt) {testIt--; ply.TestIt(); plz.TestIt(); }
01417
01418 TVectorD dW(npt21),dWy(npt),dWz(npt);
01419 TMatrixD dWW(npt21,npt21),dWWy(npt,npt),dWWz(npt,npt);
01420 double Xi2y = ply.Deriv(dWy,dWWy);
01421 double Xi2z = plz.Deriv(dWz,dWWz);
01422 double Fcn = Xi2y+Xi2z;
01423
01424
01425
01426
01427 double gi[2][3],gj[2][3],G[2][3][3];
01428 for (int ipt=0;ipt<npt;ipt++) {
01429 int idx[3]={0,ipt+1,ipt+1+npt};
01430 myDers(mDat[0],wy[ipt],wz[ipt],pt[ipt].dens,gi,G);
01431 for (int i=0;i<3;i++) {
01432 dW[idx[i]] += dWy[ipt]*gi[0][i]+dWz[ipt]*gi[1][i];
01433 for (int j=0;j<=i;j++) {
01434 dWW[idx[i]][idx[j]] += dWy[ipt]*G[0][i][j]+dWz[ipt]*G[1][i][j];
01435 } }
01436 for (int jpt=0;jpt<npt;jpt++) {
01437 int jdx[3]={0,jpt+1,jpt+1+npt};
01438 myDers(mDat[0],wy[jpt],wz[jpt],pt[jpt].dens,gj,0);
01439
01440 for (int i=0;i<3;i++) {
01441 for (int j=0;j<=i;j++) {
01442 if(idx[i]<jdx[j]) break;
01443 dWW[idx[i]][jdx[j]] += dWWy[ipt][jpt]*gi[0][i]*gj[0][j]
01444 +dWWz[ipt][jpt]*gi[1][i]*gj[1][j];
01445 } } } }
01446
01447
01448
01449 double C,Cd[3],Cdd[3][3];
01450 for (int ipt=0;ipt<npt;ipt++) {
01451 int idx[3]={0,ipt+1,ipt+1+npt};
01452 C = myFake( mDat[0],wy[ipt],wz[ipt],pt[ipt].dens,Cd,Cdd);
01453 Fcn += - log(wy[ipt]) - log(wz[ipt]) - log(C);
01454 dW[idx[1]] += -1./wy[ipt];
01455 dW[idx[2]] += -1./wz[ipt];
01456 dWW[idx[1]][idx[1]]+= 1./pow(wy[ipt],2);
01457 dWW[idx[2]][idx[2]]+= 1./pow(wz[ipt],2);
01458 for (int i=0;i<3;i++) {
01459 dW[idx[i]] += - Cd[i]/C;
01460 for (int j=0;j<=i;j++) {
01461 dWW[idx[i]][idx[j]] += (-Cdd[i][j]+Cd[i]*Cd[j]/C)/C;
01462 } } }
01463
01464
01465 for (int ii=0;ii<npt21;ii++) {
01466 if (ii && ii<=npt) wy[ii-1]*=cos2Psi[ii];
01467 dW[ii] /=cos2Psi[ii];
01468 for (int jj=0;jj<=ii;jj++) {
01469 dWW[ii][jj]/=cos2Psi[ii]*cos2Psi[jj];
01470 } }
01471
01472
01473 TVectorD wt(npt21);
01474 wt.SetSub(1,wy); wt.SetSub(1+npt,wz);
01475 for (int ii=1;ii<npt21;ii++) {
01476 double wi = wt[ii],wi2 =wi*wi;
01477 dW [ii] *= -wi2;
01478 dWW[ii][0] *= -wi2;
01479 for (int jj=1;jj<=ii;jj++) {
01480 double wj2 = wt[jj]*wt[jj];
01481 dWW[ii][jj] *= wi2*wj2;
01482 }
01483 dWW[ii][ii] += -2*dW[ii]*wi;
01484 }
01485
01486
01487
01488 for (int i=1;i<npt21;i++) {
01489 for (int j=0;j< i;j++) {
01490 assert(fabs(dWW[j][i])<=0);
01491 dWW[j][i] = dWW[i][j];
01492 } }
01493
01494 Di = toPars*dW;
01495
01496 myTCL::mxmlrtS(toPars,dWW,Dij);
01497 for (int i=0;i<mNPars;i++) {
01498 for (int j=0;j<i;j++) {
01499 assert(fabs(dWW[i][j]-dWW[j][i])<1e-10);}}
01500 return Fcn;
01501 }
01502
01503 double HitPars_t::DERIV(int npt, const MyPull *pt,TVectorD &Di,TMatrixD &Dij,int maxTrk)
01504 {
01505
01506 Di.ResizeTo(mNPars);
01507 Dij.ResizeTo(mNPars,mNPars);
01508 Di = 0.; Dij=0.; mNTrks=0;
01509 TVectorD Ti ,TiM (mNPars);
01510 TMatrixD Tij,TijM(mNPars,mNPars);
01511 int jl=0,trk=pt[0].trk;
01512 double xl=0;
01513 double sum=0,sumM=0;
01514 int NSave = (int)sqrt((npt/30.));
01515
01516 for (int jr=1; 1; jr++) {
01517 double xl0 = xl; xl = pt[jr].xyz[0];
01518 if (jr<npt && pt[jr].trk==trk && xl0<xl) continue;
01519 int n = jr-jl;
01520 if (n>15) {
01521 double tsum = Deriv(n,pt+jl,Ti,Tij);
01522 sum += tsum;
01523 Di += Ti ;
01524 Dij += Tij;
01525 mNTrks++;
01526 if (mNTrks>=maxTrk) break;
01527 if (((mNTrks+1)%NSave)==0 ) {
01528 sumM += sum; sum= 0.;
01529 TiM += Di ; Di = 0.;
01530 TijM += Dij; Dij= 0.;
01531 }
01532 }
01533 if (jr>=npt) break;
01534 trk = pt[jr].trk;
01535 jl = jr;
01536 }
01537 sum +=sumM;
01538 Di +=TiM ;
01539 Dij +=TijM;
01540 return sum;
01541 }
01542
01543 void HitPars_t::myDers( double fake,double wy, double wz, double dens
01544 , double g[2][3],double G[2][3][3])
01545 {
01546
01547
01548
01549 double u[3]={fake,wy,wz};
01550 double Cd[3],Cdd[3][3];
01551 double (*myCdd)[3] = Cdd;
01552 if (!G) myCdd=0;
01553 double C0 = myFake(fake,wy,wz,dens,Cd,myCdd);
01554
01555 memset(g[0],0,2*3*sizeof(g[0][0]));
01556 for (int iyz=0;iyz<2;iyz++) {
01557 for (int ivar=0;ivar<3;ivar++) {
01558 g[iyz][ivar] = u[iyz+1]*Cd[ivar];
01559 if (ivar==iyz+1) g[iyz][ivar]+= C0; }}
01560
01561 if (!G) return;
01562
01563 memset(G[0][0],0,2*3*3*sizeof(G[0][0][0]));
01564 for (int iyz=0;iyz<2;iyz++) {
01565 for (int i=0;i<3;i++) {
01566 for (int j=0;j<3;j++) {
01567 G[iyz][i][j] = u[iyz+1]*Cdd[i][j];
01568 if (i==iyz+1) G[iyz][i][j]+=Cd[j];
01569 if (j==iyz+1) G[iyz][i][j]+=Cd[i];
01570 } } }
01571
01572 }
01573
01574 double HitPars_t::myFake( double fake,double wy, double wz, double dens
01575 , double Cd[3],double Cdd[3][3])
01576 {
01577
01578
01579
01580 double dens2 = 2*dens;
01581
01582 double A = 3.14/sqrt(wy*wz);
01583 double Ay = -0.5*A/wy ;
01584 double Az = -0.5*A/wz ;
01585
01586 double C0 = (1+fake*dens2*A);
01587 if (!Cd) return C0;
01588 Cd[0] = dens2*A;
01589 Cd[1] = fake*dens2*Ay;
01590 Cd[2] = fake*dens2*Az;
01591 if (!Cdd) return C0;
01592
01593 double Ayy = 0.5*(-Ay + A/wy)/wy;
01594 double Ayz = 0.5*(-Az )/wy;
01595 double Azz = 0.5*(-Az + A/wz)/wz;
01596
01597 Cdd[0][0] = 0;
01598 Cdd[0][1] = dens2*Ay;
01599 Cdd[0][2] = dens2*Az;
01600 Cdd[1][0] = Cdd[0][1];
01601 Cdd[1][1] = fake*dens2*Ayy;
01602 Cdd[1][2] = fake*dens2*Ayz;
01603 Cdd[2][0] = Cdd[0][2];
01604 Cdd[2][1] = Cdd[1][2];
01605 Cdd[2][2] = fake*dens2*Azz;
01606 return C0;
01607 }
01608
01609 void HitPars_t::Print(const HitPars_t *init) const
01610 {
01611 const char *TYZ[]={"y","z"};
01612 const char *bnd;
01613
01614 printf("HitPars(Fake ) ");
01615 bnd = (Lim(0))? "*":" ";
01616 if (init) printf("\t ");
01617 printf("\tpout=%8.4g(+-%8.4g)%s\n",mDat[0],mDrr[0],bnd);
01618
01619 for (int iDet=0;iDet<mNDets;iDet++) {
01620 for (int iYZ=0;iYZ<2;iYZ++) {
01621 int ln = Len(iDet,iYZ);
01622 int ipar0 = IPar(iDet,iYZ);
01623 assert(ipar0>=0);
01624 for (int ip=0;ip<ln;ip++) {
01625 bnd = (Lim(ipar0+ip))? "*":" ";
01626 double p = mPars[iDet][iYZ][ip];
01627 double e = mErrs[iDet][iYZ][ip];
01628 printf("HitPars(%d,%s,%d) ",iDet,TYZ[iYZ],ip);
01629 if (init) {
01630 printf("\tinit=%8.4g ",init->mPars[iDet][iYZ][ip]);}
01631 printf("\tpout=%8.4g(+-%8.4g)%s",p,e,bnd);
01632 printf("\n");
01633 } } }
01634
01635 }
01636
01637 double HitPars_t::Diff(const HitPars_t &init) const
01638 {
01639 double pctMax=0;
01640 for (int i=0;i<mNPars;i++) {
01641 double dif = fabs((*this)[i]-init[i]);
01642 double err = init.Err(i);
01643 if (err<=0) err = 0.5*((*this)[i]+init[i])+1e-10;
01644 dif/= err;
01645 if (pctMax<dif) pctMax=dif;
01646 }
01647 return pctMax*100;
01648 }
01649
01650 int HitPars_t::Test(int npt, const MyPull *pt) const
01651 {
01652 enum {kNTrkTst=5};
01653 TVectorD Di ,DiT ;
01654 TMatrixD Dij,DijT;
01655 HitPars_t HP(*this);
01656 HP[0]=9;
01657 int npars = HP.NPars();
01658 double fcn = HP.DERIV(npt,pt,Di,Dij,kNTrkTst);
01659 for (int i=0;i<npars;i++) {
01660 double sav = HP[i];
01661 double delta = fabs(sav)*1e-3;
01662 if (delta<1e-6) delta=1e-6;
01663 HP[i]=sav+delta;
01664 double fcnT = HP.DERIV(npt,pt,DiT,DijT,kNTrkTst);
01665 HP[i]=sav;
01666 double ana = 0.5*(Di[i]+DiT[i]);
01667 double num = (fcnT-fcn)/delta;
01668 double dif = (num-ana)/(fabs(num+ana)+1e-5);
01669 if (fabs(dif)>0.001)
01670 printf("\nDer(%2d): ana = %g \tnum = %g \tdif = %g\n\n",i,ana,num,dif);
01671 for (int k=0;k<npars;k++) {
01672 ana = 0.5*(Dij[i][k]+DijT[i][k]);
01673 num = (DiT[k]-Di[k])/delta;
01674 dif = (num-ana)/(fabs(num+ana)+1e-5);
01675 if (fabs(dif)>0.001)
01676 printf("Der(%2d,%2d): ana = %g \tnum = %g \tdif = %g\n",i,k,ana,num,dif);
01677 }
01678
01679 }
01680 return 0;
01681 }
01682
01683 double HitPars_t::Err(const double Pars[3],int nPars,const double A[3])
01684 {
01685 static const double tenMicrons = 1e-3;
01686 static const double min2Err = tenMicrons*tenMicrons;
01687 static const double max2Err = 1.;
01688 double err = TCL::vdot(A,Pars,nPars);
01689 if (nPars<=1) return err;
01690 if (err<min2Err) err=min2Err;
01691 if (err>max2Err) err=max2Err;
01692 return err;
01693 }
01694
01695 void HitPars_t::Limit()
01696 {
01697 for (int i=0;i<mNPars;i++) {
01698 if (mDat[i]<mMin[i]) mDat[i]=mMin[i];
01699 if (mDat[i]>mMax[i]) mDat[i]=mMax[i];
01700 }
01701 }
01702 #if 1
01703
01704 void HitPars_t::Prep(int npt, const MyPull *pt,TVectorD &Y,TVectorD &Z
01705 ,TVectorD &S,TVectorD &cos2Psi)
01706 {
01707 static int myDebug=0;
01708
01709 Y.ResizeTo(npt);
01710 Z.ResizeTo(npt);
01711 S.ResizeTo(npt);
01712 if (myDebug) Show(npt,pt);
01713 int npt21 = npt*2+1;
01714 cos2Psi.ResizeTo(npt21);cos2Psi = 1.;
01715
01716 THelixFitter helx;
01717 for (int ipt=0;ipt<npt;ipt++) {
01718 TVector3 point;
01719 point[0] = pt[ipt].xhg();
01720 point[1] = pt[ipt].yhg();
01721 point[2] = pt[ipt].zhg();
01722 helx.Add(point[0],point[1],point[2]);
01723 double emx[3]={pow(pt[ipt].hye,2),0,pow(pt[ipt].hye,2)};
01724 helx.AddErr(emx,pow(pt[ipt].hze,2));
01725 }
01726 helx.Fit();
01727 THelixTrack move(helx);
01728 double s = 0;
01729 for (int ipt=0;ipt<npt;ipt++) {
01730 TVector3 point;
01731 point[0] = pt[ipt].xhg();
01732 point[1] = pt[ipt].yhg();
01733 point[2] = pt[ipt].zhg();
01734 double ds = move.Path(point[0],point[1]);
01735 if (ipt) s+=ds;
01736 move.Move(ds);
01737 TVector3 pos(move.Pos());
01738 TVector3 dir(move.Dir());
01739 TVector3 nor(dir[1],-dir[0],0.); nor = nor.Unit();
01740 double dy = (point-pos)*nor;
01741 double dz = point[2]-pos[2];
01742 S[ipt]=s;
01743 Y[ipt]=dy;
01744 Z[ipt]=dz;
01745 cos2Psi[ipt+1]=pow(cos(pt[ipt].psi),2);
01746 }
01747 }
01748 #endif//1
01749
01750 void HitPars_t::Show(int npt, const MyPull *pt)
01751 {
01752
01753
01754
01755
01756
01757
01758 static TCanvas *myCanvas=0;
01759 static TGraph *graph[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
01760
01761 float X[3][3][100],Y[3][3][100];
01762 int N[3][3];
01763
01764 if (!myCanvas) {myCanvas = new TCanvas("C1","",600,800);
01765 myCanvas->Divide(1,3);}
01766 if(pt) {
01767 double curv = 0,xPrev,yPrev; int nCurv = 0;
01768 for (int i=0;i<9;i++) {delete graph[0][i];graph[0][i]=0;}
01769 for (int ig=0;ig<3;ig++) {
01770 int n=0;
01771 double s = 0;
01772 xPrev = 3e33;
01773 for (int ipt=0;ipt<npt;ipt++) {
01774 const MyPull *node = pt+ipt;
01775
01776 double xNode = node->x_g();
01777 double yNode = node->y_g();
01778 if (xPrev<3e33) {
01779 double ds = sqrt(pow(xNode-xPrev,2)+pow(yNode-yPrev,2));
01780 double si = 0.5*ds*curv; if (si>0.99) si=0.99;
01781 if (si>0.01) ds = 2*asin(si)/curv;
01782 s += ds;
01783 }
01784 xPrev = xNode;
01785 yPrev = yNode;
01786
01787 if (ig==0) { curv += node->curv; nCurv++; continue;}
01788 if (ig==1) {
01789 X[0][ig][n] = node->x_g();
01790 Y[0][ig][n] = node->y_g();
01791 Y[2][ig][n] = node->z_g();
01792 } else {
01793 X[0][ig][n] = node->xhg();
01794 Y[0][ig][n] = node->yhg();
01795 Y[2][ig][n] = node->zhg();
01796 }
01797
01798 if (n) {
01799 float xh = X[0][ig][n]-X[0][ig][0];
01800 float yh = Y[0][ig][n]-Y[0][ig][0];
01801 float rh = xh*xh+yh*yh+1E-10;
01802 X[1][ig][n-1] = xh/rh;
01803 Y[1][ig][n-1] = yh/rh;
01804 }
01805 X[2][ig][n]=s;
01806 n++;
01807 }
01808 if (ig==0) { curv=fabs(curv)/nCurv; continue;}
01809 N[0][ig] = n;
01810 N[1][ig] = n-1;
01811 N[2][ig] = n;
01812 }
01813
01814 for (int ip=0;ip<3;ip++) {
01815 double xMin=999,xMax=-999,yMin=999,yMax=-999;
01816 for (int ig=1;ig<3;ig++) {
01817 for (int j=0;j<N[ip][ig];j++) {
01818 double x = X[ip][ig][j];
01819 if (xMin> x) xMin = x;
01820 if (xMax< x) xMax = x;
01821 double y = Y[ip][ig][j];
01822 if (yMin> y) yMin = y;
01823 if (yMax< y) yMax = y;
01824 }
01825 }
01826 X[ip][0][0] = xMin; Y[ip][0][0] = yMin;
01827 X[ip][0][1] = xMin; Y[ip][0][1] = yMax;
01828 X[ip][0][2] = xMax; Y[ip][0][2] = yMin;
01829 X[ip][0][3] = xMax; Y[ip][0][3] = yMax;
01830 N[ip][0] = 4;
01831 }
01832 static const char *opt[]={"AP","Same CP","Same *"};
01833 for (int ip=0;ip<3;ip++) {
01834 for (int ig =0;ig<3;ig++) {
01835 graph[ip][ig] = new TGraph(N[ip][ig] , X[ip][ig], Y[ip][ig]);
01836 if(ig==2) graph[ip][ig]->SetMarkerColor(kRed);
01837 myCanvas->cd(ip+1); graph[ip][ig]->Draw(opt[ig]);
01838 }
01839 }
01840
01841 }
01842
01843
01844 if (!myCanvas) return;
01845 myCanvas->Modified();
01846 myCanvas->Update();
01847 while(!gSystem->ProcessEvents()){};
01848 }
01849
01850 int HitPars_t::Test()
01851 {
01852 return 0;
01853 }
01854
01855
01856
01857
01858 double myTCL::vmaxa(const double *a,int na)
01859 {
01860 double r=0;
01861 for (int i=0;i<na;i++){if (r < fabs(a[i])) r = fabs(a[i]);}
01862 return r;
01863 }
01864
01865 double myTCL::vmaxa(const TVectorD &a)
01866 {
01867 return vmaxa(a.GetMatrixArray(),a.GetNrows());
01868 }
01869
01870 void myTCL::vfill(double *a,double f,int na)
01871 {
01872 for (int i=0;i<na;i++) {a[i]=f;}
01873 }
01874
01875 void myTCL::eigen2(const double err[3], double lam[2], double eig[2][2])
01876 {
01877
01878 double spur = err[0]+err[2];
01879 double det = err[0]*err[2]-err[1]*err[1];
01880 double dis = spur*spur-4*det;
01881 if (dis<0) dis = 0;
01882 dis = sqrt(dis);
01883 lam[0] = 0.5*(spur+dis);
01884 lam[1] = 0.5*(spur-dis);
01885 eig[0][0] = 1; eig[0][1]=0;
01886 if (dis>1e-6*spur) {
01887 if (fabs(err[0]-lam[0])>fabs(err[2]-lam[0])) {
01888 eig[0][1] = 1; eig[0][0]= -err[1]/(err[0]-lam[0]);
01889 } else {
01890 eig[0][0] = 1; eig[0][1]= -err[1]/(err[2]-lam[0]);
01891 }
01892 double tmp = sqrt(eig[0][0]*eig[0][0]+eig[0][1]*eig[0][1]);
01893 eig[0][0]/=tmp; eig[0][1]/=tmp;
01894 }
01895 eig[1][0]=-eig[0][1]; eig[1][1]= eig[0][0];
01896 }
01897
01898
01899
01900
01901
01902
01903
01904
01905
01906
01907
01908
01909
01910
01911
01912
01913
01914 double myTCL::simpson(const double *F,double A,double B,int NP)
01915 {
01916 int N2=NP-1;
01917 assert(N2>0 && !(N2&1));
01918 double S1=F[N2-1];
01919 double S2=0;
01920
01921 for (int N = 1;N<=N2-3;N+=2) {S1+=F[N];S2+=F[N+1];}
01922 S1=S1+S1+S2;
01923 double H=(F[0]+F[N2]+S1+S1)*(B-A)/(3*N2);
01924 return H;
01925 }
01926
01927 void myTCL::mxmlrt(const TMatrixD &A,const TMatrixD &B,TMatrixD &X)
01928 {
01929 int nRowA = A.GetNrows();
01930 int nColA = A.GetNcols();
01931 int nRowB = B.GetNrows();
01932 int nColB = B.GetNcols(); if(nColB){}
01933 assert(nColA ==nRowB);
01934 X.ResizeTo(nRowA,nRowA);
01935 TCL::mxmlrt(A.GetMatrixArray(),B.GetMatrixArray()
01936 ,X.GetMatrixArray(),nRowA,nColA);
01937
01938 }
01939
01940 void myTCL::mxmlrtS(const double *A,const double *B,double *X,int nra,int nca)
01941 {
01942 TCL::vzero(X,nra*nra);
01943 for (int i=0,ii=0;i<nra;i++,ii+=nca) {
01944 for (int j=0,jj=0;j<nca;j++,jj+=nca) {
01945 if(!A[ii+j]) continue;
01946 for (int k=0,kk=0;k<=i;k++,kk+=nca) {
01947 double &Xik =X[i*nra+k];
01948 for (int l=0;l<nca;l++) {
01949 if(!A[kk+l]) continue;
01950 Xik +=A[ii+j]*A[kk+l]*B[jj+l];
01951 } } } }
01952 for (int i=0;i<nra;i++){
01953 for (int k=0;k<i ;k++){X[k*nra+i] = X[i*nra+k];}}
01954 }
01955
01956 void myTCL::mxmlrtS(const TMatrixD &A,const TMatrixD &B,TMatrixD &X)
01957 {
01958 int nRowA = A.GetNrows();
01959 int nColA = A.GetNcols();
01960 int nRowB = B.GetNrows();
01961 int nColB = B.GetNcols(); if(nColB){}
01962 assert(nColA ==nRowB);
01963 X.ResizeTo(nRowA,nRowA);
01964 myTCL::mxmlrtS(A.GetMatrixArray(),B.GetMatrixArray()
01965 ,X.GetMatrixArray(),nRowA,nColA);
01966
01967 }
01968
01969 double myTCL::vasum(const double *a, int na)
01970 {
01971 double sum = 0;
01972 for (int i=0;i<na;i++) { sum += fabs(a[i]);}
01973 return sum;
01974 }
01975
01976 double myTCL::vasum(const TVectorD &a)
01977 {
01978 return vasum(a.GetMatrixArray(),a.GetNrows());
01979 }
01980
01981 int myTCL::SqProgSimple( TVectorD &x
01982 ,const TVectorD &g,const TMatrixD &G
01983 ,const TVectorD &Min
01984 ,const TVectorD &Max, int iAktp)
01985 {
01986 static int nCall=0; nCall++;
01987 enum {kINIT=0,kADDCONS,kFREECONS};
01988 int kase = kINIT;
01989 int nPars = g.GetNrows();
01990 TVectorD xx(x),gg(g),add(nPars);
01991 TArrayI Side(nPars);
01992 int nCons=0,addCons = -1,freCons=-1,freSide=0,addSide=0,con=0;
01993 double maxGra=3e33;
01994 int iAkt = iAktp;
01995 while(1946) {
01996
01997 freCons=-1; freSide=0;
01998 if (nCons && kase==kFREECONS ) {
01999 double tryGra=kSMALL; freCons=-1;
02000 for (int ix=0;ix<nPars;ix++) {
02001 if(!Side[ix]) continue;
02002 double gra = gg[ix]*Side[ix];
02003 if (gra< tryGra) continue;
02004 if (gra>=maxGra) continue;
02005 freCons=ix; tryGra=gra;
02006 }
02007 if (freCons>=0) {
02008 maxGra = tryGra;
02009 freSide = Side[freCons];
02010 Side[freCons]=0;
02011 nCons--;
02012 }
02013 }
02014
02015 if(kase==kFREECONS && freCons<0) { break;}
02016
02017
02018 TMatrixD S(G);
02019 TVectorD B(gg);
02020 if (nCons) {
02021 for (int ix=0;ix<nPars;ix++) {
02022 if (Side[ix]==0) continue;
02023 for (int jx=0;jx<nPars;jx++) {S[ix][jx]=0; S[jx][ix]=0;}
02024 S[ix][ix]=1; B[ix]=0;
02025 } }
02026 if (iAkt==0 ) {
02027
02028 double det=S.Determinant();
02029 if (fabs(det)<1e-100) return -99;
02030 S.Invert(0);
02031
02032
02033 add = (-1.)*(S*B);
02034 double along = B*add;
02035 if (along>0) add*=(-1.);
02036 }
02037
02038 if (iAkt==1 ) {
02039 double bb = (B*B);
02040 double bSb = (B*(S*B));
02041 double tau = -bb/(fabs(bSb)+3e-33);
02042 add = tau*B;
02043
02044 }
02045 if(kase==kFREECONS && freSide) {
02046 if (add[freCons]*freSide > -kSMALL) {
02047 Side[freCons]=freSide; nCons++; continue;}
02048 }
02049
02050
02051 double fak=1;
02052 addCons = -1; addSide = 0;
02053 con = 0;
02054 for (int ix=0;ix<nPars;ix++) {
02055 if (Side[ix]) {add[ix]=0; con = 100*con+ix+1;continue;}
02056 double xi = xx[ix]+fak*add[ix];
02057 if (xi < Min[ix]){fak = (Min[ix]-xx[ix])/add[ix]; addCons=ix;addSide=-1;}
02058 if (xi > Max[ix]){fak = (Max[ix]-xx[ix])/add[ix]; addCons=ix;addSide= 1;}
02059 assert(fak<=1. && fak>=0.);
02060 }
02061 add*=fak;
02062 xx+= add;
02063 gg += G*add;
02064 maxGra=3e33;
02065 kase = kFREECONS; if (!addSide) continue;
02066 kase = kADDCONS;
02067 xx[addCons] = (addSide<0)? Min[addCons]:Max[addCons];
02068
02069 Side[addCons] = addSide ;nCons++;
02070
02071 }
02072 x = xx;
02073 return abs(con);
02074 }
02075 #endif
