StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
fiterr.C
1 //#ifndef __CINT__
2 #include <stdio.h>
3 #include <iostream>
4 #include <fstream>
5 #include "TSystem.h"
6 #include "TMath.h"
7 #include "TBenchmark.h"
8 #include "TFile.h"
9 #include "TTree.h"
10 #include "TH1.h"
11 #include "TF1.h"
12 #include "TCanvas.h"
13 #include "TGraph.h"
14 #include "TCernLib.h"
15 #include "TMatrixD.h"
16 #include "TRandom.h"
17 #include "TVectorD.h"
18 #include "TVector2.h"
19 #include "TVector3.h"
20 #include "StarRoot/TTreeIter.h"
21 #include "TTable.h"
22 #include "TInterpreter.h"
23 #include "TPolinom.h"
24 #include "THelixTrack.h"
25 #include <vector>
26 //#endif
27 #define SQ(X) ((X)*(X))
28 #define PW2(X) ((X)*(X))
29 #define PW3(X) ((X)*(X)*(X))
30 
31 enum {kNOBAGR=0,kNDETS=4,kMINHITS=50};
32 static const double WINDOW_NSTD=3;
33 
34 static const double kAGREE=1e-7,kSMALL=1e-9,kBIG=1.,kDAMPR=0.1;
35 static const int MinErr[4][2] = {{200,200},{200,200},{30,500},{20,20}};
36 static char gTimeStamp[16]={0};
37 class MyPull;
38 class HitPars_t;
39 int fiterr(const char *opt);
40 double Fit(HitPars_t &pout);
41 void AveRes();
42 double AveErr();
43 int kind(double xl);
44 void DbInit();
45 void DbDflt();
46 void DbEnd();
47 void Init(HitPars_t &hitp);
48 void HInit();
49 void HFit();
50 void HEnd();
51 void CheckErr();
52 void FillPulls(int befAft);
53 double newPull(int iYZ,const MyPull& myRes,const HitPars_t &pars);
54 TBranch *tb=0;
55 TFile *pulFile =0;
56 TTree *pulTree=0;
57 TH1F *hz=0;
58 TH1F *hx=0;
59 TH1F *hh[100];
60 TCanvas *C[10];
61 
62 class myTCL {
63 public:
64 static double vmaxa (const double *a, int na);
65 static double vmaxa (const TVectorD &a);
66 static void vfill (double *a,double f,int na);
67 static void mxmlrt(const TMatrixD &A,const TMatrixD &B,TMatrixD &X);
68 static void mxmlrtS(const TMatrixD &A,const TMatrixD &B,TMatrixD &X);
69 static void mxmlrtS(const double *A,const double *B,double *X,int nra,int nca);
70 static void eigen2(const double err[3], double lam[2], double eig[2][2]);
71 static double simpson(const double *F,double A,double B,int NP);
72 static double vasum(const double *a, int na);
73 static double vasum(const TVectorD &a);
74 static int SqProgSimple( TVectorD &x
75  ,const TVectorD &g,const TMatrixD &G
76  ,const TVectorD &Min
77  ,const TVectorD &Max,int iAkt);
78 };
79 
80 
81 class HitAccr {
82 public:
83 int iDet;
84 double A[2][3];
85 double Pull[2];
86 double PredErr[2];
87 };
88 
89 
90 class MyPull {
91 public:
92 float x_g() const { return grf*cos(gpf);}
93 float y_g() const { return grf*sin(gpf);}
94 float z_g() const { return gzf ;}
95 
96 float xhg() const { return xyz[0]*cos(ang) - xyz[1]*sin(ang);}
97 float yhg() const { return xyz[0]*sin(ang) + xyz[1]*cos(ang);}
98 float zhg() const { return xyz[2] ;}
99 public:
100 int trk;
101 float xyz[3];
102 float psi;
103 float dip;
104 float pt;
105 float ypul;
106 float zpul;
107 float yfit;
108 float zfit;
109 float pye; //pull y errvmaxa
110 float pze; //pull z err
111 float uyy; //untouched y err
112 float uzz; //untouched z err
113 float fye; //fited y err
114 float fze; //fited z err
115 float hye; //hited y err
116 float hze; //hited z err
117 float curv; //curva
118 float dens;
119 float grf; // Rxy of Fit in global Sti frame
120 float gpf; // Phi of Fit in global Sti frame
121 float gzf; // Z of Fit in global Sti frame
122 float ang; // rotation angle
123 };
124 
125 class HitPars_t {
126 public:
127 HitPars_t();
128 HitPars_t(const HitPars_t &fr);
129 const double &Err(int idx) const {return mDrr[idx];}
130  double &Err(int idx) {return mDrr[idx];}
131 double Err( int iYZ,const HitAccr &accr) const;
132 double Err(int iDet,int iYZ,const double A[3]) const;
133 const double &operator[](int i) const;
134  double &operator[](int i);
135 HitPars_t &operator= (const HitPars_t &fr);
136 HitPars_t &operator* (double f);
137 HitPars_t &operator+=(const double *f);
138 HitPars_t &operator= (double f);
139 HitPars_t &operator+=(const TVectorD &add);
140 int NPars() const {return mNPars ;}
141 int Len(int iDet,int iYZ=0) const {return mLen[iDet][iYZ];}
142 int Lim(int i) const ;
143 const double &Min(int i) const {return mMin[i];}
144  double &Min(int i) {return mMin[i];}
145 const double &Max(int i) const {return mMax[i];}
146  double &Max(int i) {return mMax[i];}
147 
148 const double &Min(int iDet,int iYZ,int i) const {return mMin[IPar(iDet,iYZ)+i];}
149 const double &Max(int iDet,int iYZ,int i) const {return mMax[IPar(iDet,iYZ)+i];}
150  double &Min(int iDet,int iYZ,int i) {return mMin[IPar(iDet,iYZ)+i];}
151  double &Max(int iDet,int iYZ,int i) {return mMax[IPar(iDet,iYZ)+i];}
152 
153 int IPar(int iDet,int iYZ,int *npars=0) const;
154 void Set(int iDet,int iYZ,int nini,const double *ini);
155 void Limit();
156 double Deriv(int npt, const MyPull *pt,TVectorD &Di,TMatrixD &Dij) const;
157 double DERIV(int npt, const MyPull *pt,TVectorD &Di,TMatrixD &Dij,int maxTrk=9999999);
158 int Test(int npt, const MyPull *pt) const;
159 void Print(const HitPars_t *init=0) const;
160 double Diff (const HitPars_t &init) const;
161 static void Prep(int npt, const MyPull *pt,TVectorD &Y,TVectorD &Z
162  ,TVectorD &S,TVectorD &cos2Psi);
163 static int Test();
164 static void HitCond(const MyPull& myRes,HitAccr &acc);
165 static void Show(int npt,const MyPull *pt);
166 static double Dens(double rxy,int ntk);
167 static double Err(const double Pars[3],int nPars,const double A[3]);
168 static void myDers( double fake, double wy, double wz, double dens
169  , double g[2][3],double G[2][3][3]);
170 static double myFake( double fake, double wy, double wz, double dens
171  , double Cd[3],double Cdd[3][3]);
172 public:
173 int mNDets;
174 int mNPars;
175 int mNTrks;
176 int mLen[kNDETS][2];
177 double *mPars[kNDETS][2]; //mPars[iDet][iYZ][3]
178 double *mErrs[kNDETS][2]; //mErrs[iDet][iYZ][3]
179 double mDat[1+kNDETS*2*3];
180 double mDrr[1+kNDETS*2*3];
181 double mMin[1+kNDETS*2*3];
182 double mMax[1+kNDETS*2*3];
183 
184 private:
185 static void myTrans( double fake,double wy, double wz, double dens, TMatrixD *mx);
186 };
187 HitPars_t operator+(const HitPars_t &a,const double *add);
188 HitPars_t operator+(const HitPars_t &a,const TVectorD &add);
189 
190 
191 
192 
193 //______________________________________________________________________________
194 class Poli2 {
195 public:
196 enum {kXMAX=200};
197 Poli2(int npw=1);
198 Poli2(int npw,int n,const double *X,const double *Y,const double *W);
199 void Init();
200 void Clear();
201 void Add(double x, double y, double w);
202 double Xi2() const {return fXi2;}
203 double Xi2(int ipt) const;
204 double EvalXi2() const;
205 int Pw() const {return fPw ;}
206 int NPt() const {return fN ;}
207 double Fun( double x ) const;
208 double dXi2dW( int ipt ) const;
209 double d2Xi2d2W( int ipt,int jpt ) const;
210 double Deriv( TVectorD &Di, TMatrixD &Dij ) const;
211 void TestIt() const;
212 void MyTest() const;
213 static void Test();
214 private:
215 int fPw;
216 char fBeg[1];
217 int fN;
218 double fX[kXMAX];
219 double fY[kXMAX];
220 double fW[kXMAX];
221 double fXi2;
222 double fX0,fY0;
223 char fEnd[1];
224 TVectorD fP;
225 TVectorD fB;
226 TMatrixD fA;
227 TMatrixD fAi;
228 };
229 
230 
231 
232 #if !defined(__MAKECINT__)
233 
234 
235 
236 
237 
238 std::vector<MyPull> MyVect;
239  double aveRes[4][6],aveTrk[4][2][3];
240  int numRes[4];
241  double pSTI[8][3] = {{0.000421985 ,0.00124327 ,0.0257111 }
242  ,{0.000402954 ,0.00346896 ,0.0377259 }
243  ,{0. ,0.0009366 ,0.0004967 }
244  ,{0.00018648 ,0.00507244 ,0.002654 }
245  ,{0.0030*0.0030 ,0.0 ,0.0 }
246  ,{0.0700*0.0700 ,0.0 ,0.0 }
247  ,{0.0080*0.0080 ,0.0 ,0.0 }
248  ,{0.0080*0.0080 ,0.0 ,0.0 }};
249 HitPars_t HitErr;
250 
251 static const char *DETS[]={"OutY","OutZ","InnY","InnZ"
252  ,"SsdY","SsdZ","SvtY","SvtZ"};
253 static const char *DETZ[]={"Out" ,"Inn","Ssd","Svt"};
254  int FitOk[4]={0,0,0,0};
255 static char dbFile[4][100] = {
256 "StarDb/Calibrations/tracker/tpcOuterHitError.20050101.235959.C",
257 "StarDb/Calibrations/tracker/tpcInnerHitError.20050101.235959.C",
258 "StarDb/Calibrations/tracker/ssdHitError.20050101.235959.C" ,
259 "StarDb/Calibrations/tracker/svtHitError.20050101.235959.C" };
260 
261 static TTable *dbTab[4];
262 void myBreak(int kase) {
263  static int myKase=-1946;
264  if (kase != myKase) return;
265  printf("BOT OHO\n");
266 }
267 
268 int fiterr(const char *opt)
269 {
270  if (!opt) opt = "H";
271  int optH = strstr(opt,"H")!=0;
272  int optU = strstr(opt,"U")!=0;
273  int optT = strstr(opt,"T")!=0;
274  memcpy(gTimeStamp,strstr(opt,"20"),15);
275 
276  DbInit();
277  if (optH) HInit();
278  TTreeIter th(pulTree);
279  th.AddFile("pulls.root");
280 // th.ls();
281 // return;
282  const int &run = th("mRun");
283  const int &evt = th("mEvt");
284 //const int *&mNTrks = th("mNTrks[2]");
285  const int &mTrksG = th("mTrksG");
286  const int &nGlobs = th("mHitsG");
287  const float *&mCurv = th("mHitsG.mCurv");
288  const float *&mPt = th("mHitsG.mPt");
289  const short *&mTrackNumber = th("mHitsG.mTrackNumber");
290  const float *&mNormalRefAngle = th("mHitsG.mNormalRefAngle");
291  const UChar_t *&nTpcHits = th("mHitsG.nTpcHits");
292  const UChar_t *&nSsdHits = th("mHitsG.nSsdHits");
293  const UChar_t *&nSvtHits = th("mHitsG.nSvtHits");
294  const float *&lYPul = th("mHitsG.lYPul");
295  const float *&lZPul = th("mHitsG.lZPul");
296  const float *&lXHit = th("mHitsG.lXHit");
297  const float *&lYHit = th("mHitsG.lYHit");
298  const float *&lZHit = th("mHitsG.lZHit");
299  const float *&lYFit = th("mHitsG.lYFit");
300  const float *&lZFit = th("mHitsG.lZFit");
301  const float *&lPsi = th("mHitsG.lPsi");
302  const float *&lDip = th("mHitsG.lDip");
303  const float *&lYFitErr = th("mHitsG.lYFitErr");
304  const float *&lZFitErr = th("mHitsG.lZFitErr");
305  const float *&lYHitErr = th("mHitsG.lYHitErr");
306  const float *&lZHitErr = th("mHitsG.lZHitErr");
307  const float *&lYPulErr = th("mHitsG.lYPulErr");
308  const float *&lZPulErr = th("mHitsG.lZPulErr");
309  const float *&gRFit = th("mHitsG.gRFit");
310  const float *&gPFit = th("mHitsG.gPFit");
311  const float *&gZFit = th("mHitsG.gZFit");
312 
313  printf("*lYPul = %p\n",lYPul);
314  printf("&run=%p &evt=%p\n",&run,&evt);
315  MyPull mp;
316  MyVect.clear();
317  int nTpc=0,nPre=0,iTk=-1,iTkSkip=-1;
318  while (th.Next())
319  {
320 // printf("serial = %d run=%d evt=%d nGlo = %d\n",n,run,evt,nGlobs);n++;
321  float rxy=0;
322  for (int ih=0;ih<nGlobs;ih++) {
323  if (mTrackNumber[ih]==iTkSkip) continue;
324  float rxy00 = rxy; rxy=lXHit[ih];
325  if (mTrackNumber[ih]!=iTk || rxy<rxy00) {
326  iTk = mTrackNumber[ih];
327  iTkSkip = iTk;
328  if (nTpcHits[ih] <25) continue;
329  if (fabs(mCurv[ih])>1./200) continue;
330  if (fabs(mPt[ih]/cos(lDip[ih]))<0.5) continue;
331  int nSS = nSsdHits[ih]+nSvtHits[ih];
332  if (nSS && nSS<2) continue;
333  if (!nSS && nPre && nTpc>100 && nTpc>3*nPre) continue;
334  nTpc += (nSS==0);
335  nPre += (nSS!=0);
336  iTkSkip = -1;
337  }
338 
339  memset(&mp,0,sizeof(mp));
340  mp.trk = mTrackNumber[ih];
341  mp.pye = lYPulErr[ih];
342  mp.pze = lZPulErr[ih];
343  float uyy = (lYPulErr[ih]-lYHitErr[ih])*(lYHitErr[ih]+lYPulErr[ih]);
344 // printf ("yy = %g %g %g \n",yy,lYHitErr[ih],lYPulErr[ih]);
345 // if (uyy <1e-8 ) continue;
346  float uzz = (lZPulErr[ih]-lZHitErr[ih])*(lZHitErr[ih]+lZPulErr[ih]);
347 // printf ("zz = %g\n",zz);
348 // if (uzz <1e-8 ) continue;
349  mp.uyy = uyy;
350  mp.uzz = uzz;
351  mp.fye = lYFitErr[ih];
352  mp.fze = lZFitErr[ih];
353  mp.hye = lYHitErr[ih];
354  mp.hze = lZHitErr[ih];
355  mp.xyz[0] = lXHit[ih];
356  if (mp.xyz[0]<4) continue;
357  mp.xyz[1] = lYHit[ih];
358  mp.xyz[2] = lZHit[ih];
359  mp.psi = lPsi[ih] ;
360  mp.dip = lDip[ih] ;
361  mp.ypul = lYPul[ih];
362 // if (fabs(mp.ypul)<1e-8) continue;
363  mp.zpul = lZPul[ih];
364 // if (fabs(mp.zpul)<1e-8) continue;
365  mp.yfit = lYFit[ih];
366  mp.zfit = lZFit[ih];
367  mp.dens = HitPars_t::Dens(mp.xyz[0],mTrksG);
368  mp.grf = gRFit[ih];
369  mp.gpf = gPFit[ih];
370  mp.gzf = gZFit[ih];
371  mp.curv = mCurv[ih];
372  mp.ang = mNormalRefAngle[ih];
373 
374  MyVect.push_back(mp);
375 
376  }
377 // printf("lYPul[0]=%g\n",lYPul[0]);
378  }
379  printf("fiterr:: %d hits accepted\n\n",MyVect.size());
380  printf("fiterr:: %d Tpc and %d Svt/Ssd tracks accepted\n\n",nTpc,nPre);
381  AveRes();
382  HFit();
383  DbDflt();
384  Init(HitErr);
385  if (optH) CheckErr();
386  if (optH) FillPulls(0);
387  if (optT) return HitErr.Test(MyVect.size(),&(MyVect[0]));
388 
389  double maxPct = Fit(HitErr);
390  maxPct= AveErr();
391  if (optH) FillPulls(1);
392  if (optU) DbEnd();
393  if (optH) HEnd();
394 //return (maxPct<10) ? 99:0;
395  return (maxPct< 3) ? 99:0;
396 }
397 //______________________________________________________________________________
398 Poli2::Poli2(int npw):fP(npw+1),fB(npw+1),fA(npw+1,npw+1),fAi(npw+1,npw+1)
399 {
400 assert(0);
401  fPw = npw;
402  Clear();
403 
404 }
405 //______________________________________________________________________________
406 void Poli2::Clear()
407 {
408  memset(fBeg,0,fEnd-fBeg+1);
409 }
410 //______________________________________________________________________________
411 Poli2::Poli2(int npw,int N,const double *X,const double *Y,const double *W)
412 :fP(npw+1),fB(npw+1),fA(npw+1,npw+1),fAi(npw+1,npw+1)
413 {
414 assert(0);
415  fPw = npw;
416  Clear();
417  fN = N;
418  int nb = N*sizeof(double);
419  memcpy(fX,X,nb);
420  memcpy(fY,Y,nb);
421  memcpy(fW,W,nb);
422 }
423 //______________________________________________________________________________
424 void Poli2::Add(double x, double y, double w)
425 {
426  fX[fN]=x; fY[fN]=y, fW[fN]=w; fN++;
427  assert(fN<=kXMAX);
428 }
429 //______________________________________________________________________________
430 void Poli2::Init()
431 {
432  double Wt=0,Wx=0,Wy=0;
433  for (int i=0;i<fN;i++) { Wt += fW[i]; Wx += fW[i]*fX[i];Wy += fW[i]*fY[i];}
434  fX0 = Wx/Wt; fY0 = Wy/Wt;
435  for (int i=0;i<fN;i++) { fX[i]-=fX0; fY[i]-=fY0;}
436 
437  double A[3][3]={{0}},B[3]={0};
438  double x[3]={1};
439  fXi2=0;
440  for (int i=0;i<fN;i++) {
441  double w = fW[i];
442  x[1] = fX[i]; x[2]=x[1]*x[1];
443  double y = fY[i],y2=y*y;
444  for (int j=0;j<=fPw;j++) { B[j] +=w*x[j]*y;
445  for (int k=0;k<=j ;k++) { A[j][k] +=w*x[j]*x[k];}}
446  fXi2 +=w*y2;
447  }
448  for (int j=0;j<=fPw;j++){ fB[j] = B[j];
449  for (int k=0;k<=j ;k++){ fA(j,k) = A[j][k]; fA(k,j)=A[j][k];}}
450 
451  double det=0;
452  fAi=fA;fAi.InvertFast(&det);
453  fP=fAi*fB;
454  fXi2 -= (fB*fP);
455 
456 }
457 //______________________________________________________________________________
458 double Poli2::Fun( double x ) const
459 {
460  x -= fX0;
461  double xx[3]={1,x,x*x};
462  TVectorD XX(fPw+1,xx);
463  return fP*XX + fY0;
464 }
465 //______________________________________________________________________________
466 double Poli2::Xi2(int ipt) const
467 {
468  double dy = Fun(fX[ipt]+fX0) - (fY[ipt]+fY0);
469  return dy*dy*fW[ipt];
470 }
471 //______________________________________________________________________________
472 double Poli2::EvalXi2() const
473 {
474  double sum=0;
475  for (int ipt=0;ipt<fN;ipt++) {sum+=Xi2(ipt);}
476  return sum;
477 }
478 //______________________________________________________________________________
479 double Poli2::dXi2dW( int ipt ) const
480 {
481 // Xi2 = B*Ai*B
482 // dXi2 = 2*dB*Ai*B - B*Ai*dA*Ai*B
483 // dXi2 = 2*dB*Ai*B - P*dA*P
484 
485  double der = fY[ipt]*fY[ipt];
486  TVectorD dB(fPw+1);
487  TMatrixD dA(fPw+1,fPw+1);
488  double x[3]={1,fX[ipt],fX[ipt]*fX[ipt]};
489  for (int j=0;j<=fPw;j++) { dB[j] = x[j]*fY[ipt];
490  for (int k=0;k<=fPw;k++) { dA(j,k) = x[j]*x[k]; }}
491  der -= (2*(dB*(fAi*fB))-(fP*(dA*fP)));
492  return der;
493 }
494 //______________________________________________________________________________
495 //______________________________________________________________________________
496 double Poli2::d2Xi2d2W( int ipt,int jpt ) const
497 {
498 // P = Ai*B
499 // dXi2/dWi = 2*diB*Ai*B - B*Ai*diA*Ai*B
500 // d2Xi2/dWi/dWj = -2*diB*Ai*djA*Ai*B
501 // +2*diB*Ai*djB
502 // -2*djB*Ai*diA*Ai*B
503 // +2*B*Ai*djA*Ai*diA*Ai*B
504 // d2Xi2/dWi/dWj = -2*diB*Ai*djA*P
505 // +2*diB*Ai*djB
506 // -2*djB*Ai*diA*P
507 // +2*P*djA*Ai*diA*P
508 
509  double xi[3]={1,fX[ipt],fX[ipt]*fX[ipt]};
510  double xj[3]={1,fX[jpt],fX[jpt]*fX[jpt]};
511 
512  TVectorD diB(fPw+1) ,djB(fPw+1);
513  TMatrixD diA(fPw+1,fPw+1),djA(fPw+1,fPw+1);
514  for (int j=0;j<=fPw;j++) { diB[j]=fY[ipt]*xi[j]; djB[j]=fY[jpt]*xj[j];
515  for (int k=0;k<=fPw;k++) { diA(j,k) = xi[j]*xi[k]; djA(j,k) = xj[j]*xj[k];}}
516  double der = -2*(-((fAi*diB)*(djA*fP))
517  +(diB*(fAi*djB))
518  -((fAi*djB)*(diA*fP))
519  +(fP*(djA*(fAi*(diA*fP)))));
520  return der;
521 }
522 //______________________________________________________________________________
523 double Poli2::Deriv( TVectorD &Di, TMatrixD &Dij ) const
524 {
525  Di.ResizeTo(fN);
526  Dij.ResizeTo(fN,fN);
527  TVectorD diB(fPw+1) ,djB(fPw+1);
528  TMatrixD diA(fPw+1,fPw+1),djA(fPw+1,fPw+1);
529  for (int ipt=0;ipt<fN;ipt++) {
530  double der = fY[ipt]*fY[ipt];
531  double xi[3]={1,fX[ipt],fX[ipt]*fX[ipt]};
532  for (int k=0;k<=fPw;k++) { diB[k] = fY[ipt]*xi[k];
533  for (int l=0;l<=fPw;l++) { diA(k,l) = xi[k]*xi[l];}}
534  der -= (2*(diB*(fAi*fB))-(fP*(diA*fP)));
535  Di[ipt] = der;
536  for (int jpt=0;jpt<=ipt;jpt++) {
537  double xj[3]={1,fX[jpt],fX[jpt]*fX[jpt]};
538  for (int k=0;k<=fPw;k++) { djB[k] = fY[jpt]*xj[k];
539  for (int l=0;l<=fPw;l++) { djA(k,l) = xj[k]*xj[l];}}
540  der = -2*(-((fAi*diB)*(djA*fP))
541  +(diB*(fAi*djB))
542  -((fAi*djB)*(diA*fP))
543  +(fP*(djA*(fAi*(diA*fP)))));
544  Dij[ipt][jpt]=der;
545  Dij[jpt][ipt]=der;
546  }
547  }
548  return Xi2();
549 }
550 
551 //______________________________________________________________________________
552 void Poli2::MyTest() const
553 {
554  printf("Print B\n");
555  fB.Print();
556  printf("Print A\n");
557  fA.Print();
558  printf("Print Ai\n");
559  fAi.Print();
560  printf("Print A*Ai\n");
561  TMatrixD tst(fPw+1,fPw+1);
562  tst = fA*fAi;
563  tst.Print();
564 }
565 //______________________________________________________________________________
566 void Poli2::Test()
567 {
568  int npw=2;
569  double X[20],Y[20],W[20];
570  for (int i=0;i<20;i++) {
571  X[i]=i;
572  Y[i]= 3+X[i]*(.02+.03*X[i]);
573  W[i]= 1-0.01*i;
574  Y[i]+=gRandom->Gaus(0,sqrt(1./W[i]));
575  }
576  Poli2 pp(npw,20,X,Y,W);
577  pp.Init();
578 // pp.MyTest();
579  double Xi2 = pp.Xi2();
580  double Xi2Eva = pp.Xi2();
581  printf ("Xi2 = %g == %g\n",Xi2,Xi2Eva);
582  for (int i=0; i<20; i++) {
583  printf(" %g %g : %g\n",X[i],Y[i],pp.Fun(X[i]));
584  }
585 // check d/dW
586  printf("\n check d/dWi\n");
587  double delta = 1e-3;
588  for (int k=0;k<20;k++) {
589  Poli2 ppp(npw);
590  for (int i=0;i<20;i++) {
591  double w = W[i];
592  if (i==k) w+=delta;
593  ppp.Add(X[i],Y[i],w);
594  }
595  ppp.Init();
596  double ana = pp.dXi2dW(k);
597  double num = (ppp.Xi2()-pp.Xi2())/delta;
598  double dif = 2*(num-ana)/(fabs(num+ana)+3e-33);
599  if (fabs(dif)<0.01) continue;
600  printf ("dXi2dW(%2d) \tana=%g \tnum = %g \tdif=%g\n",k,ana,num,dif);
601  }
602 // check d2/dWi/DWj
603  printf("\n check d2/dWi/DWj\n");
604  for (int i=0;i<20;i++) {
605  double sav = W[i]; W[i]=sav+delta;
606  Poli2 ppp(npw,20,X,Y,W);
607  W[i]=sav;
608  ppp.Init();
609  for (int j=0;j<20;j++) {
610  double ana = pp.d2Xi2d2W(i,j);
611  double num = (ppp.dXi2dW(j)-pp.dXi2dW(j))/delta;
612  double dif = 2*(num-ana)/(fabs(num+ana)+3e-33);
613  if (fabs(dif)<0.01) continue;
614  printf ("d2Xi2d2W(%2d,%2d) \tana=%g \tnum = %g \tdif=%g\n",i,j,ana,num,dif);
615  } }
616  TVectorD g;
617  TMatrixD gg;
618  pp.Deriv(g,gg);
619  for (int i=0;i<20;i++) {
620  double tst = pp.dXi2dW(i)-g[i];
621  if (fabs(tst)>1e-10) printf("g[%d] %g **************\n",i,tst);
622  for (int j=0;j<20;j++) {
623  double tst = pp.d2Xi2d2W(i,j)-gg[i][j];
624  if (fabs(tst)>1e-10) printf("gg[%d][%d] %g **************\n",i,j,tst);
625  } }
626 
627 
628 }
629 //______________________________________________________________________________
630 void Poli2::TestIt() const
631 {
632  TVectorD g,gT;
633  TMatrixD G,GT;
634  Poli2 pp(*this);
635  pp.Init();
636  double Xi2 = pp.Deriv(g,G);
637 
638  double delta = 1e-3;
639  for (int k=0;k<fN;k++) {
640  double sav = pp.fW[k];
641  pp.fW[k]=sav+delta;
642  pp.Init();
643  double Xi2T = pp.Deriv(gT,GT);
644  pp.fW[k]=sav;
645 
646  double ana = 0.5*(g[k]+gT[k]);
647  double num = (Xi2T-Xi2)/delta;
648  double dif = 2*(num-ana)/(fabs(num+ana)+1e-10);
649  if (fabs(dif)>0.01)
650  printf ("\ndXi2dW(%2d) \tana=%g \tnum = %g \tdif=%g\n\n",k,ana,num,dif);
651 // check d2/dWi/DWj
652  for (int i=0;i<fN;i++) {
653  ana = 0.5*(GT[k][i]+G[k][i]);
654  num = (gT[i]-g[i])/delta;
655  dif = 2*(num-ana)/(fabs(num+ana)+1e-10);
656  if (fabs(dif)>0.01)
657  printf ("d2Xi2dW2(%2d,%2d) \tana=%g \tnum = %g \tdif=%g\n",k,i,ana,num,dif);
658  } }
659 
660 
661 }
662 //______________________________________________________________________________
663 void poli2()
664 {
665 Poli2::Test();
666 }
667 //______________________________________________________________________________
668 void Init(HitPars_t &hiterr)
669 {
670  hiterr[0] = 0.0;
671  for (int iDet=0;iDet<kNDETS;iDet++) {
672  if (numRes[iDet]<kMINHITS) continue;
673  int n = (iDet<=1)? 3:1;
674  hiterr.Set(iDet,0,n,pSTI[iDet*2+0]);
675  hiterr.Set(iDet,1,n,pSTI[iDet*2+1]);
676  hiterr.Min(iDet,0,0) = pow(MinErr[iDet][0]*1.e-4,2);
677  hiterr.Min(iDet,1,0) = pow(MinErr[iDet][1]*1.e-4,2);
678  }
679 }
680 //______________________________________________________________________________
681 //______________________________________________________________________________
682 class FitState_t : public HitPars_t {
683  public:
684  FitState_t();
685  FitState_t(HitPars_t &hp);
686  const HitPars_t &Pars() const { return (const HitPars_t &)(*this);}
687  HitPars_t &Pars() { return ( HitPars_t &)(*this);}
688  int Saddle() const {return neg<0;}
689  int LimX(int i) const;
690  double Der() const {return der;}
691  double Fcn() const {return fcn;}
692  void Deriv(const std::vector<MyPull> &MyVect);
693  int MaxStp(TVectorD &add,int mode) const;
694  void MakeErrs();
695  int operator<(const FitState_t &other) const;
696 FitState_t &operator=(const FitState_t &other);
697 static int FixWeak( TVectorD &g, TMatrixD &G);
698  public:
699  TVectorD g;
700  TMatrixD G;
701  char myBeg[1];
702  int npt;
703  int ider;
704  double fak;
705  double fcn;
706  double der;
707  double neg;
708  char myEnd[1];
709 };
710 //______________________________________________________________________________
711 FitState_t::FitState_t()
712 {
713  memset(myBeg,0,myEnd-myBeg+1);
714  der = 1e99;fcn = 1e99;neg =-1e99;
715 }
716 //______________________________________________________________________________
717 FitState_t::FitState_t(HitPars_t &hp) :HitPars_t(hp)
718 {
719  memset(myBeg,0,myEnd-myBeg+1);
720  der = 1e99;fcn = 1e99;neg =-1e99;
721 }
722 //______________________________________________________________________________
723 void FitState_t::Deriv(const std::vector<MyPull> &MyVect)
724 {
725  npt = MyVect.size();
726  fcn = DERIV(npt,&(MyVect[0]),g,G);
727  fcn/=npt;
728  der = -1; ider =-1; neg = 0;
729  for (int i=0;i<mNPars;i++) {
730  if (G[i][i]<neg) neg =G[i][i];
731  if(LimX(i)) continue;
732  if (fabs(g[i])>der) { ider = i; der = fabs(g[i]);}
733  }
734  der/=npt;
735  double myMax=myTCL::vmaxa(G.GetMatrixArray(),G.GetNoElements());
736  fak = pow(2.,-int(log(myMax)/log(2.)));
737  G*=fak; g*=fak;
738  if (neg>=0) { //Check positiveness
739  for (int i = 0; i<mNPars && neg>0; i++){
740  for (int j = 0; j<i && neg>0; j++) {
741  if (pow(G[i][j],2) >= G[i][i]*G[j][j]) neg = -1;
742  } }
743  }
744 }
745 
746 
747 //______________________________________________________________________________
748 int FitState_t::operator<(const FitState_t &old) const
749 {
750 //if (Pos()<old.Pos()) return 1;
751  double delta = (Fcn()-old.Fcn())/Fcn();
752 
753  if (Der()<0 ) return 0;
754  if (Fcn()<old.Fcn()) return 1;
755  if (Der()<old.Der()
756  && Fcn()-old.Fcn()-0.1*fabs(old.Fcn())<0) return 1;
757  return 0;
758 }
759 //______________________________________________________________________________
760 FitState_t &FitState_t::operator=(const FitState_t &other)
761 {
762  memcpy(myBeg,other.myBeg,myEnd-myBeg+1);
763  g.ResizeTo(other.g);
764  G.ResizeTo(other.G);
765  g = other.g;
766  G = other.G;
767  HitPars_t::operator=(other);
768  return *this;
769 }
770 //______________________________________________________________________________
771 int FitState_t::MaxStp(TVectorD &add,int mode) const
772 {
773  int lim=0;
774  double fak=1;
775  for (int i=0;i<mNPars;i++) {
776  double maxStp = (0.01+mDat[i])*0.3;
777  if (maxStp>fak*fabs(add[i])) continue;
778  lim = lim*100+i+1;
779  if (!mode) { add[i] = (add[i]<0)? -maxStp:maxStp;}
780  else { fak = maxStp/fabs(add[i]);}
781  }
782  if (fak<1) add*=fak;
783  return lim;
784 }
785 //______________________________________________________________________________
786 void FitState_t::MakeErrs()
787 {
788 // Evaluate errors
789  TMatrixD E(G);
790  for (int i=0;i<mNPars;i++) {
791  if (!Lim(i)) continue;
792  for (int j=0;j<mNPars;j++) {E(i,j) = 0; E(j,i) = 0;}; E(i,i) = 1;
793  }
794  double det=12345; E.InvertFast(&det);
795  E*=fak;
796  for (int i=0;i<mNPars;i++) {Err(i) = sqrt(E(i,i));}
797 }
798 //______________________________________________________________________________
799 int FitState_t::LimX(int i) const
800 {
801  int lim = Lim(i);
802  if (!lim) return 0;
803  if (lim<0 && g[i]<0) return 0;
804  if (lim>0 && g[i]>0) return 0;
805  return lim;
806 }
807 //______________________________________________________________________________
808 int FitState_t::FixWeak( TVectorD &g, TMatrixD &G)
809 {
810  int n = g.GetNrows();
811  int nfix = 0;
812  double ave = myTCL::vasum(g)/n;
813  for (int i=0;i<n;i++) {
814  if (fabs(g[i])>=ave) continue;
815  nfix++; g[i]=0;
816  for (int j=0;j<n;j++) {G[i][j]=0; G[j][i]=0;}
817  G[i][i]=1;
818  }
819  return nfix;
820 }
821 //______________________________________________________________________________
822 double Fit(HitPars_t &pout)
823 {
824  enum {kMAXITER=10000};
825 static int const MAXCUT[2]={4,10};
826 static int kount=0;
827 
828  int ifail = 99,nPars,iAkt=0,iter,nCut,lim=0,con=0;
829  nPars = pout.NPars();
830  HitPars_t init(pout);
831  init.Limit();
832 
833  TVectorD add(nPars),g(nPars),ge(nPars),Gg(nPars);
834  FitState_t Best(init),Now(init);
835  double dif,difFcn=0,difDer;
836 
837  static int idebug = 1;
838 
839  iAkt=0;nCut=0;
840 
841  for (iter=0;iter<kMAXITER;iter++) {// Iterations
842  kount++;
843  Now.Deriv(MyVect);
844  difFcn = Now.fcn-Best.fcn;
845  difDer = Now.der-Best.der;
846  if (idebug) {
847  printf("Fit(%3d) \tFcn = %g(%g) \tDer(%2d)=%g(%g) \tLim=%d \tCon=%d\tAkt=%d Cut=%d Sad=%d\n"
848  ,iter,Now.fcn,difFcn,Now.ider,Now.der,difDer,lim,con,iAkt,nCut,Now.Saddle());
849  }
850  if (Now.Der()<0) {ifail=1; break;
851  } else if (Now.Der() < kAGREE) {
852  Best=Now; ifail=0; break;
853 
854  } else if (Now < Best ) {//Accept
855  nCut = 0; Best = Now; iAkt = 0;
856 
857  } else if (nCut>MAXCUT[iAkt]) {//Cut step failed
858  if (iAkt==0) Now=Best;
859  Best=Now; nCut=0; iAkt=1-iAkt;
860 
861  } else {
862  nCut++;
863  add*= 0.5;
864  Now.Pars() = Best.Pars()+add;
865  continue;
866  }
867 
868 // SOLVING
869 
870 // Sometimes try grad method randomly
871 static int rnd=0; rnd++;
872  if (!(rnd&3)) iAkt=1;
873 
874 // if (Now.Saddle() && iter&1)iAkt=1;
875  TVectorD P0(nPars,&Best[0]);
876  for (int jAkt=iAkt;jAkt<2;jAkt++) {
877  iAkt = jAkt;
878  TVectorD P1(P0);
879  con = myTCL::SqProgSimple(P1,Now.g,Now.G
880  ,TVectorD(nPars,&Now.Min(0))
881  ,TVectorD(nPars,&Now.Max(0)),jAkt);
882  if (con<0) continue;
883  add = P1-P0;
884  double along = -(add*Now.g);
885  if (along <0) continue;
886  lim=Best.MaxStp(add,1);
887  Now.Pars() = Best.Pars()+add;
888  break;
889  }
890  }// End Iters
891 
892  if (ifail==0 || ifail==99) Best.MakeErrs();
893  pout=Best.Pars();
894 
895  dif = pout.Diff(init);
896  printf("\nFit: Iter=%d Fcn=%g Der=%g Dif=%6.3g%% Fail=%d\n"
897  ,iter,Best.fcn,Best.der,dif,ifail);
898  pout.Print(&init);
899  return (ifail)? 1e10:dif;
900 
901 }
902 
903 //______________________________________________________________________________
904 void AveRes()
905 {
906  memset(aveRes[0] ,0,sizeof(aveRes));
907  memset(aveTrk[0][0],0,sizeof(aveTrk));
908  memset(numRes ,0,sizeof(numRes));
909  for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
910  MyPull &myRes = MyVect[jhit];
911  int jdx = kind(myRes.xyz[0]);
912  aveRes[jdx][0]+= myRes.ypul*myRes.ypul;
913  aveRes[jdx][1]+= myRes.zpul*myRes.zpul;
914  aveRes[jdx][2]+= myRes.uyy;
915  aveRes[jdx][3]+= myRes.uzz;
916  aveRes[jdx][4]+= pow(myRes.xyz[1]-myRes.yfit,2);
917  aveRes[jdx][5]+= pow(myRes.xyz[2]-myRes.zfit,2);
918  if (hh[jdx*2+0+30]) hh[jdx*2+0+30]->Fill(myRes.ypul);
919  if (hh[jdx*2+1+30]) hh[jdx*2+1+30]->Fill(myRes.zpul);
920  numRes[jdx]++;
921  HitAccr ha;
922  HitPars_t::HitCond(myRes,ha);
923  TCL::vadd(aveTrk[jdx][0],ha.A[0],aveTrk[jdx][0],3);
924  TCL::vadd(aveTrk[jdx][1],ha.A[1],aveTrk[jdx][1],3);
925 
926  }
927 // CleanUp
928  int ihit=0;
929  for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
930  MyPull &myRes = MyVect[jhit];
931  int jdx = kind(myRes.xyz[0]);
932  if (numRes[jdx]<kMINHITS) continue;
933  aveRes[jdx][0]+= myRes.ypul*myRes.ypul;
934  if (ihit!=jhit) MyVect[ihit]=MyVect[jhit];
935  ihit++;
936  }
937  MyVect.resize(ihit);
938 
939  for(int jdx=0;jdx<4;jdx++) {
940  if (numRes[jdx]<kMINHITS) continue;
941  double f = 1./numRes[jdx];
942  TCL::vscale(aveRes[jdx],f,aveRes[jdx],6);
943  TCL::vscale(aveTrk[jdx][0],f,aveTrk[jdx][0],6);
944  double hitYErr,hitZErr;
945  TCL::vsub(aveRes[jdx]+0,aveRes[jdx]+2,aveRes[jdx]+2,2);
946  hitYErr = aveRes[jdx][2];
947  hitYErr = (hitYErr<0)? -sqrt(-hitYErr):sqrt(hitYErr);
948  hitZErr = aveRes[jdx][3];
949  hitZErr = (hitZErr<0)? -sqrt(-hitZErr):sqrt(hitZErr);
950  printf("AveRes:: N=%5d \taveRes[%s][yz]=%g %g \tavePul[yz]=%g %g \thitErr(yz)=%g %g\n\n"
951  ,numRes[jdx],DETZ[jdx]
952  ,sqrt(aveRes[jdx][4]),sqrt(aveRes[jdx][5])
953  ,sqrt(aveRes[jdx][0]),sqrt(aveRes[jdx][1])
954  ,hitYErr,hitZErr);
955 
956 
957 
958  }
959 }
960 //______________________________________________________________________________
961 int kind(double x)
962 {
963 static const double radds[5]={120,25,16,0,0};
964  for(int jdx=0;1;jdx++) {if (x>radds[jdx]) return jdx;}
965 }
966 //___________________________________myRes___________________________________________
967 void HInit()
968 {
969  memset(hh,0,sizeof(hh));
970  memset(C,0,sizeof(C));
971 // Pulls before and after
972 static const char *NamPuls[]={
973 "YOut.bef","YOut.aft","ZOut.bef","ZOut.aft",
974 "YInn.bef","YInn.aft","ZInn.bef","ZInn.aft",
975 "YSsdBef","YSsdAft","ZSsdBef","ZSsdAft",
976 "YSvtBef","YSvtAft","ZSvtBef","ZSvtAft"};
977 static const char *NamRes[]={
978 "YOut.res","ZOut.res",
979 "YInn.res","ZInn.res",
980 "YSsd.res","ZSsd.res",
981 "YSvt.res","ZSvt.res"};
982 
983  for (int ih=0;ih<16;ih++) {
984  hh[ih] = new TH1F(NamPuls[ih],NamPuls[ih],100,-3,3);}
985  C[0] = new TCanvas("C0","",600,800);
986  C[0]->Divide(1,8);
987  C[1] = new TCanvas("C1","",600,800);
988  C[1]->Divide(1,8);
989  for (int ih=0;ih<8;ih++) {C[0]->cd(ih+1);hh[ih+0]->Draw();}
990  for (int ih=0;ih<8;ih++) {C[1]->cd(ih+1);hh[ih+8]->Draw();}
991 
992  hh[20]= new TH1F("ChekErr","ChekErr",100,-100.,100.);
993  hh[21]= new TH1F("ChekPul","ChekPul",100,-100.,100.);
994  C[2] = new TCanvas("C2","",600,800);
995  C[2]->Divide(1,2);
996  for (int ih=0;ih<2;ih++) {C[2]->cd(ih+1);hh[ih+20]->Draw();}
997 // hh[20]->Draw();
998 
999  for (int ih=0;ih<8;ih++) {
1000  hh[ih+30] = new TH1F(NamRes[ih],NamRes[ih],100,-0.3,0.3);}
1001  C[3] = new TCanvas("C3","",600,800);
1002  C[3]->Divide(1,8);
1003  for (int ih=0;ih<8;ih++) {C[3]->cd(ih+1);hh[ih+30]->Draw();}
1004 
1005 }
1006 
1007 //______________________________________________________________________________
1008 void FillPulls(int befAft)
1009 {
1010  MyPull myRes;
1011  for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
1012  myRes = MyVect[jhit];
1013  int jdx = kind(myRes.xyz[0]);
1014  if (befAft==0) {
1015  hh[jdx*4+0]->Fill(myRes.ypul/(myRes.pye));
1016  hh[jdx*4+2]->Fill(myRes.zpul/(myRes.pze));
1017  } else {
1018  hh[jdx*4+1]->Fill(newPull(0,myRes,HitErr));
1019  hh[jdx*4+3]->Fill(newPull(1,myRes,HitErr));
1020  }
1021  }
1022 }
1023 //______________________________________________________________________________
1024 double newPull(int iYZ,const MyPull& myRes,const HitPars_t &pars)
1025 {
1026  HitAccr accr;
1027  HitPars_t::HitCond(myRes,accr);
1028  double S = pars.Err(iYZ,accr)+accr.PredErr[iYZ];
1029  if (S<1e-6) S=1e-6;
1030  return accr.Pull[iYZ]/sqrt(S);
1031 }
1032 
1033 //______________________________________________________________________________
1034 void DbInit()
1035 {
1036  gSystem->Load("libStDb_Tables.so");
1037  TString command;
1038  TTable *newdat = 0;
1039  for (int idb=0;idb<4;idb++) {
1040  memcpy(strstr(dbFile[idb],"20"),gTimeStamp,15);
1041  int ready=0;
1042  if (!gSystem->AccessPathName(dbFile[idb])) {//file exists
1043  command = ".L "; command += dbFile[idb];
1044  gInterpreter->ProcessLine(command);
1045  newdat = (TTable *) gInterpreter->Calc("CreateTable()");
1046  command.ReplaceAll(".L ",".U ");
1047  gInterpreter->ProcessLine(command);
1048  ready = 2006;
1049  } else {
1050  newdat = (TTable *)gInterpreter->Calc("new St_HitError(\"someHitError\",1)");
1051  newdat->SetUsedRows(1);
1052  }
1053  assert(newdat);
1054  dbTab[idb] = newdat;
1055  }
1056 }
1057 //______________________________________________________________________________
1058 void DbDflt()
1059 {
1060 // Set initial values if db file was non existing
1061 
1062  for (int idb=0;idb<4;idb++)
1063  {
1064  if (numRes[idb]<kMINHITS) continue;
1065  double *d = (double *)dbTab[idb]->GetArray();
1066  if (d[0]<=0) {d[0]=aveRes[idb][0];d[3]=aveRes[idb][1];}
1067  TCL::ucopy(d+0,pSTI[idb*2+0],3);
1068  TCL::ucopy(d+3,pSTI[idb*2+1],3);
1069  }
1070 }
1071 //______________________________________________________________________________
1072 void DbEnd()
1073 {
1074  double par[2][3];
1075  for (int idb=0;idb<kNDETS;idb++)
1076  {
1077  memset(par,0,sizeof(par));
1078  if (!HitErr.mPars[idb][0]) continue;
1079  if (!HitErr.Len(idb)) continue;
1080  TString ts(dbFile[idb]);
1081  ts +=".BAK";
1082  gSystem->Rename(dbFile[idb],ts);
1083  int ny = HitErr.Len(idb,0);
1084  TCL::ucopy(HitErr.mPars[idb][0],par[0],ny);
1085  int nz = HitErr.Len(idb,1);
1086  TCL::ucopy(HitErr.mPars[idb][1],par[1],nz);
1087 
1088  TCL::ucopy(par[0],(double*)dbTab[idb]->GetArray()+0,3+3);
1089 
1090  std::ofstream ofs(dbFile[idb]);
1091  dbTab[idb]->SavePrimitive(ofs);
1092  }
1093 }
1094 
1095 //______________________________________________________________________________
1096 void HEnd()
1097 {
1098  for (int ic=0;ic<(int)(sizeof(C)/sizeof(void*));ic++) {
1099  TCanvas *cc = C[ic]; if (!cc) continue;
1100  cc->Modified();cc->Update();
1101  }
1102  while(!gSystem->ProcessEvents()){};
1103 }
1104 
1105 //______________________________________________________________________________
1106 void CheckErr()
1107 {
1108  for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
1109  MyPull &myRes = MyVect[jhit];
1110  int jdx = kind(myRes.xyz[0]); if(jdx){};
1111  if (jdx != 2) continue;
1112  float pye = myRes.pye;
1113  float hye = myRes.hye;
1114 
1115  float ypul = myRes.ypul;
1116  if (fabs(ypul/pye)<2) {
1117  double tmp = (ypul*ypul-pye*pye)/(hye*hye);
1118  hh[20]->Fill(tmp);
1119  }
1120  hh[21]->Fill(pye/hye);
1121  }
1122 }
1123 //______________________________________________________________________________
1124 double AveErr()
1125 {
1126 // idx=0 outer resY
1127 // idx=1 outer resZ
1128 // idx=2 inner resY
1129 // idx=3 inner resZ
1130 // idx=4 Ssd resY
1131 // idx=5 Ssd resZ
1132 // idx=4 Svt resY
1133 // idx=5 Svt resZ
1134  double pctMax=0;
1135  for (int iDet=0;iDet<kNDETS;iDet++) {
1136  for (int iYZ=0;iYZ<2;iYZ++) {
1137  int idx = iYZ + 2*iDet;
1138  int nPars = HitErr.Len(iDet,iYZ);
1139  double sum=0,sum2=0,sig[3]={0};
1140  int nTot=0;
1141  for (int jhit=0; jhit<(int)MyVect.size(); jhit++) {
1142  MyPull &myRes = MyVect[jhit];
1143  int jdx = kind(myRes.xyz[0]);
1144  if (jdx != iDet) continue;
1145  nTot++;
1146  HitAccr accr;
1147  HitPars_t::HitCond(myRes,accr);
1148  double S = HitErr.Err(iYZ,accr);
1149  double delta = pow((&myRes.ypul)[iYZ],2);
1150  double P0 = delta/((&myRes.uyy)[iYZ]+pow((&myRes.hye)[iYZ],2));
1151  double P1 = delta/((&myRes.uyy)[iYZ]+HitPars_t::Err(pSTI[idx],nPars,accr.A[iYZ]));
1152  double P2 = delta/((&myRes.uyy)[iYZ]+S);
1153  sum +=S; sum2 += S*S; sig[0]+=P0; sig[1]+=P1;sig[2]+=P2;
1154  }
1155  if (!nTot) continue;
1156  sum/=nTot; sum2/=nTot; sum2-=sum*sum; if (sum2<0) sum2=0; sum2=sqrt(sum2);
1157  sum = sqrt(sum); sum2 /= 2*sum;
1158  for (int i=0;i<3;i++) {sig[i]=sqrt(sig[i]/nTot);}
1159 
1160  int isig = int(10000*sum +0.5);
1161  int esig = int(10000*sum2+0.5);
1162  double newErr = sum;
1163  double oldErr = sqrt(HitPars_t::Err(pSTI[idx],nPars,aveTrk[iDet][iYZ]));
1164  double pct = 200*fabs(newErr-oldErr)/(newErr+oldErr);
1165  if (pct>pctMax) pctMax=pct;
1166  int iold = int(10000*oldErr +0.5);
1167 // printf("AveErr(%s): <Err.old> = %4d <Err.new> %4d(+-%4dmu) <Sig> =%g %g %g\t//evts=%d\n"
1168 // ,DETS[idx],iold,isig,esig,sig[0],sig[1],sig[2],nTot);
1169  printf("AveErr(%s): <Err.old> = %4d <Err.new> %4d(+-%4dmu) Dif=%4.1g%%\t// evts=%d\n"
1170  ,DETS[idx],iold,isig,esig,pct,nTot);
1171  } }
1172  printf("AveErr(): maxPct = %4.1g%%\n\n", pctMax );
1173 
1174 
1175  return pctMax;
1176 }
1177 //______________________________________________________________________________
1178 void HFit()
1179 {
1180  if (!hh[30+4]) return;
1181 }
1182 //_____________________________________________________________________________
1183 //_____________________________________________________________________________
1184 //_____________________________________________________________________________
1185 //_____________________________________________________________________________
1186 HitPars_t::HitPars_t ()
1187 {
1188  memset(this,0,sizeof(HitPars_t));
1189  mNPars=1;
1190  int n = sizeof(mMin)/sizeof(*mMin);
1191  myTCL::vfill(mMin,kSMALL,n);
1192  myTCL::vfill(mMax,kBIG ,n);
1193  mMax[0] = 3;
1194  mDat[0] = 0.1;
1195 }
1196 //_____________________________________________________________________________
1197 HitPars_t::HitPars_t (const HitPars_t &fr)
1198 {
1199  *this = fr;
1200 }
1201 //_____________________________________________________________________________
1202 HitPars_t &HitPars_t::operator=(const HitPars_t &fr)
1203 {
1204  memcpy(this,&fr,sizeof(HitPars_t));
1205  int inc = (char*)mDat-(char*)fr.mDat;
1206  int n = sizeof(mPars)/sizeof(mPars[0][0]);
1207  char **p = (char**)mPars[0];
1208  for (int i=0;i<n;i++) {if (p[i]) p[i] += inc;}
1209  assert(mDat <=mPars[0][0]);
1210  return *this;
1211 }
1212 //_____________________________________________________________________________
1213 const double &HitPars_t::operator[](int i) const
1214 { return mDat[i];}
1215 
1216 //_____________________________________________________________________________
1217  double &HitPars_t::operator[](int i)
1218 { return mDat[i];}
1219 //_____________________________________________________________________________
1220 HitPars_t &HitPars_t::operator*(double f)
1221 { for (int i=0;i<mNPars;i++) { (*this)[i]*=f;} return *this; }
1222 //_____________________________________________________________________________
1223 HitPars_t &HitPars_t::operator+=(const double *f)
1224 { for (int i=0;i<mNPars;i++) {(*this)[i]+=f[i];} return *this;}
1225 //_____________________________________________________________________________
1226 HitPars_t &HitPars_t::operator=(double f)
1227 { for (int i=0;i<mNPars;i++) {(*this)[i]=f;} return *this;}
1228 //_____________________________________________________________________________
1229 HitPars_t operator+(const HitPars_t &a,const double *add)
1230 {
1231  HitPars_t tmp(a); tmp+=add; return tmp;
1232 }
1233 //_____________________________________________________________________________
1234 HitPars_t operator+(const HitPars_t &a,const TVectorD &add)
1235 {
1236  HitPars_t tmp(a); tmp+=add; return tmp;
1237 }
1238 //_____________________________________________________________________________
1239 int HitPars_t::IPar(int iDet,int iYZ, int *npars) const
1240 {
1241  if (npars) *npars=Len(iDet,iYZ);
1242  int ans = mPars[iDet][iYZ]-mDat;
1243  assert(ans>=0 && ans < 100);
1244  return ans;
1245 }
1246 //_____________________________________________________________________________
1247 int HitPars_t::Lim(int i) const
1248 {
1249  int lim = 0;
1250  if (mDat[i] <= mMin[i]*1.1) lim = -1;
1251  if (mDat[i] >= mMax[i]*0.9) lim = 1;
1252  return lim;
1253 }
1254 //_____________________________________________________________________________
1255 double HitPars_t::Err(int iDet,int iYZ,const double A[3]) const
1256 {
1257  return Err(mPars[iDet][iYZ],Len(iDet,iYZ),A);
1258 }
1259 //_____________________________________________________________________________
1260 double HitPars_t::Err(int iYZ,const HitAccr &accr) const
1261 {
1262  return Err(accr.iDet,iYZ,accr.A[iYZ]);
1263 }
1264 //_____________________________________________________________________________
1265 HitPars_t &HitPars_t::operator+=(const TVectorD &add)
1266 {
1267  for (int i=0;i<mNPars;i++) {
1268  mDat[i]+=add[i];
1269  if (mDat[i]<mMin[i]) mDat[i]=mMin[i];
1270  if (mDat[i]>mMax[i]) mDat[i]=mMax[i];
1271  }
1272  return *this;
1273 }
1274 //_____________________________________________________________________________
1275 void HitPars_t::Set(int iDet,int iYZ,int nini,const double *ini)
1276 {
1277  mLen[iDet][iYZ]=nini;
1278  if (iDet+1>mNDets) mNDets=iDet+1;;
1279  mPars[iDet][iYZ]=mDat+mNPars;
1280  mErrs[iDet][iYZ]=mDrr+mNPars;
1281  mNPars+=nini;
1282  for (int i=0;i<nini;i++) {
1283  int ipar = IPar(iDet,iYZ);
1284  assert(ipar>=0);
1285  mDat[ipar+i]=ini[i];
1286  if (mDat[ipar+i]< mMin[ipar+i]) mDat[ipar+i]= mMin[ipar+i];
1287  if (mDat[ipar+i]> mMax[ipar+i]) mDat[ipar+i]= mMax[ipar+i];
1288  }
1289 }
1290 
1291 //_____________________________________________________________________________
1292 double HitPars_t::Dens(double rxy,int ntk)
1293 {
1294 // pseudo-rapidity = -log(tan(theta/2)) theta angle btw track & beam
1295 
1296 
1297  return ntk/(4*3.14*rxy*rxy);
1298 }
1299 //______________________________________________________________________________
1300 void HitPars_t::HitCond(const MyPull& myRes,HitAccr &acc)
1301 {
1302  memset(acc.A[0],0,sizeof(acc.A));
1303  acc.iDet = kind(myRes.xyz[0]);
1304  for (int iyz=0;iyz<2;iyz++) {
1305  acc.A[iyz][0]=1;
1306  acc.A[iyz][1]= 0.01*(200-fabs(myRes.xyz[2]));
1307  double ang;
1308  if (!iyz) {
1309  ang = myRes.psi; acc.Pull[0]=myRes.ypul;
1310  acc.PredErr[0] = myRes.uyy;
1311  } else {
1312  ang = myRes.dip; acc.Pull[1]=myRes.zpul;
1313  acc.PredErr[1] = myRes.uzz;
1314  }
1315  double ca=cos(ang);
1316  if (ca<0.01) ca=0.01;
1317  double ca2 = ca*ca;
1318  double ta2=(1-ca2)/ca2;
1319  acc.A[iyz][1]/= ca2;
1320  acc.A[iyz][2] = ta2;
1321  }
1322 }
1323 //______________________________________________________________________________
1324 double HitPars_t::Deriv(int npt, const MyPull *pt,TVectorD &Di,TMatrixD &Dij) const
1325 {
1326 // Ai = ErrYi**2; Bi = ErrZi**2;
1327 // Wy_i = 1./Ai ; Wz_i = 1./Bi ;
1328 // Area_i = Pi*sqrt(Ai*Bi);
1329 // -log(Fcn_i) = 0.5*(1+2*fake*area_i*dens_i)*(Wy_i*(Yi-Pol_y(Li))**2
1330 // +Wz_i*(Zi-Pol_z(Li))**2)
1331 // - 0.5*log(1+2*fake*area_i*dens_i)
1332 // - 0.5*log(Wyi) - 0.5*log(Wzi)
1333 // Fcn = Fcn_0+Fcn_1+...
1334 
1335 
1336  Di.ResizeTo(mNPars);
1337  Dij.ResizeTo(mNPars,mNPars);
1338  int npt21 = npt*2+1;
1339 
1340  HitAccr ha;
1341  Poli2 ply(2);
1342  Poli2 plz(1);
1343  TVectorD wy(npt),wz(npt),cos2Psi(npt21);
1344  TMatrixD toPars(mNPars,npt21);
1345  toPars[0][0]=1;
1346  cos2Psi=1.;
1347  TPoliFitter pfy(2),pfz(1); //?????
1348  TCircleFitter cf;
1349 #define NEWPOLIFIT
1350 #ifndef NEWPOLIFIT
1351  double len = 0,oldCurv,oldXyz[3],nowXyz[3];
1352  assert(pt[0].xyz[0]<pt[npt-1].xyz[0]);
1353  for (int ipt=0;ipt<npt;ipt++) {
1354  HitCond(pt[ipt],ha);
1355  int iDet = ha.iDet;
1356  if (!mPars[iDet][0]) continue;
1357  double nowCurv=pt[ipt].curv;
1358  nowXyz[0] = pt[ipt].grf*cos(pt[ipt].gpf);
1359  nowXyz[1] = pt[ipt].grf*sin(pt[ipt].gpf);
1360  if (ipt) { //calc length
1361  double dl = sqrt(SQ(nowXyz[0]-oldXyz[0])+SQ(nowXyz[1]-oldXyz[1]));
1362  double curv=0.5*fabs(nowCurv+oldCurv);
1363  if (dl*curv<1e-3) dl*=(1.+SQ(dl*curv)/24);
1364  else dl = 2*asin(0.5*dl*curv)/curv;
1365  len+=dl;
1366  }
1367 
1368  double dy = pt[ipt].xyz[1]-pt[ipt].yfit;
1369  double cosPsi = cos(pt[ipt].psi);
1370  cos2Psi[ipt+1] = cosPsi*cosPsi;
1371  dy *= cosPsi;
1372  double dz = pt[ipt].xyz[2]-pt[ipt].zfit;
1373  wy[ipt] = (1./Err(0,ha))/cos2Psi[ipt+1];
1374  wz[ipt] = (1./Err(1,ha));
1375  double fake = myFake( mDat[0],wy[ipt],wz[ipt],pt[ipt].dens,0,0);
1376  ply.Add(len,dy,wy[ipt]*fake);
1377  plz.Add(len,dz,wz[ipt]*fake);
1378 // pfy.Add(len,dy,1./(wy[ipt]*fake));//????
1379 // pfz.Add(len,dz,1./(wz[ipt]*fake));//????
1380 // double emx[3]={1./(wy[ipt]*fake),0,1./(wy[ipt]*fake)};
1381 // cf.Add(len,dy,emx);//????
1382  int n,ipar;
1383  for (int jk=0;jk<2;jk++) {
1384  int ip = 1 +ipt + jk*npt;
1385  ipar = IPar(iDet,jk,&n);
1386  assert(ipar>=0);
1387  for (int in=0;in<n;in++) {toPars[ipar+in][ip]+=ha.A[jk][in];}
1388  }
1389  memcpy(oldXyz,nowXyz,sizeof(oldXyz));
1390  oldCurv=nowCurv;
1391  }
1392 #endif //0
1393 #ifdef NEWPOLIFIT
1394  TVectorD Y,Z,S;
1395  Prep(npt,pt,Y,Z,S,cos2Psi);
1396  for (int ipt=0;ipt<npt;ipt++) {
1397  HitCond(pt[ipt],ha);
1398  wy[ipt] = (1./Err(0,ha))/cos2Psi[ipt+1];
1399  wz[ipt] = (1./Err(1,ha));
1400  double fake = myFake( mDat[0],wy[ipt],wz[ipt],pt[ipt].dens,0,0);
1401  ply.Add(S[ipt],Y[ipt],wy[ipt]*fake);
1402  plz.Add(S[ipt],Z[ipt],wz[ipt]*fake);
1403  pfy.Add(S[ipt],Y[ipt],1./(wy[ipt]*fake));
1404  pfz.Add(S[ipt],Z[ipt],1./(wz[ipt]*fake));
1405  double emx[3]={1./(wy[ipt]*fake),0,1./(wy[ipt]*fake)};
1406  cf.Add(S[ipt],Y[ipt],emx);
1407  int n,ipar,iDet = ha.iDet;
1408 static int myCall=0; myCall++;
1409  for (int jk=0;jk<2;jk++) {
1410  int ip = 1 +ipt + jk*npt;
1411  ipar = IPar(iDet,jk,&n);
1412  assert(ipar>=0);
1413  for (int in=0;in<n;in++) {toPars[ipar+in][ip]+=ha.A[jk][in];}
1414  }
1415  }
1416 #endif //0
1417  ply.Init();
1418  plz.Init();
1419 
1420 static int testIt=1;
1421  if (testIt) {testIt--; ply.TestIt(); plz.TestIt(); }
1422 
1423  TVectorD dW(npt21),dWy(npt),dWz(npt);
1424  TMatrixD dWW(npt21,npt21),dWWy(npt,npt),dWWz(npt,npt);
1425  double Xi2y = ply.Deriv(dWy,dWWy);
1426  double Xi2z = plz.Deriv(dWz,dWWz);
1427  double Fcn = Xi2y+Xi2z;
1428 
1429 // account that W was multiplied by fake, where
1430 // fake = (1+2*mFake*pt[ipt].dens*Pi/sqrt(wx[ipt]*wy[ipt]))
1431 
1432  double gi[2][3],gj[2][3],G[2][3][3];
1433  for (int ipt=0;ipt<npt;ipt++) {
1434  int idx[3]={0,ipt+1,ipt+1+npt};
1435  myDers(mDat[0],wy[ipt],wz[ipt],pt[ipt].dens,gi,G);
1436  for (int i=0;i<3;i++) {
1437  dW[idx[i]] += dWy[ipt]*gi[0][i]+dWz[ipt]*gi[1][i];
1438  for (int j=0;j<=i;j++) {
1439  dWW[idx[i]][idx[j]] += dWy[ipt]*G[0][i][j]+dWz[ipt]*G[1][i][j];
1440  } }
1441  for (int jpt=0;jpt<npt;jpt++) {
1442  int jdx[3]={0,jpt+1,jpt+1+npt};
1443  myDers(mDat[0],wy[jpt],wz[jpt],pt[jpt].dens,gj,0);
1444 
1445  for (int i=0;i<3;i++) {
1446  for (int j=0;j<=i;j++) {
1447  if(idx[i]<jdx[j]) break;
1448  dWW[idx[i]][jdx[j]] += dWWy[ipt][jpt]*gi[0][i]*gj[0][j]
1449  +dWWz[ipt][jpt]*gi[1][i]*gj[1][j];
1450  } } } }
1451 
1452 
1453 // account log() terms
1454  double C,Cd[3],Cdd[3][3];
1455  for (int ipt=0;ipt<npt;ipt++) {
1456  int idx[3]={0,ipt+1,ipt+1+npt};
1457  C = myFake( mDat[0],wy[ipt],wz[ipt],pt[ipt].dens,Cd,Cdd);
1458  Fcn += - log(wy[ipt]) - log(wz[ipt]) - log(C);
1459  dW[idx[1]] += -1./wy[ipt];
1460  dW[idx[2]] += -1./wz[ipt];
1461  dWW[idx[1]][idx[1]]+= 1./pow(wy[ipt],2);
1462  dWW[idx[2]][idx[2]]+= 1./pow(wz[ipt],2);
1463  for (int i=0;i<3;i++) {
1464  dW[idx[i]] += - Cd[i]/C;
1465  for (int j=0;j<=i;j++) {
1466  dWW[idx[i]][idx[j]] += (-Cdd[i][j]+Cd[i]*Cd[j]/C)/C;
1467  } } }
1468 
1469 // account Wy = wy/cos2Psi
1470  for (int ii=0;ii<npt21;ii++) {
1471  if (ii && ii<=npt) wy[ii-1]*=cos2Psi[ii];
1472  dW[ii] /=cos2Psi[ii];
1473  for (int jj=0;jj<=ii;jj++) {
1474  dWW[ii][jj]/=cos2Psi[ii]*cos2Psi[jj];
1475  } }
1476 
1477 // change W to E where W=1/E
1478  TVectorD wt(npt21);
1479  wt.SetSub(1,wy); wt.SetSub(1+npt,wz);
1480  for (int ii=1;ii<npt21;ii++) {
1481  double wi = wt[ii],wi2 =wi*wi;
1482  dW [ii] *= -wi2;
1483  dWW[ii][0] *= -wi2;
1484  for (int jj=1;jj<=ii;jj++) {
1485  double wj2 = wt[jj]*wt[jj];
1486  dWW[ii][jj] *= wi2*wj2;
1487  }
1488  dWW[ii][ii] += -2*dW[ii]*wi;
1489  }
1490 
1491 
1492 // Symmetrisation
1493  for (int i=1;i<npt21;i++) {
1494  for (int j=0;j< i;j++) {
1495  assert(fabs(dWW[j][i])<=0);
1496  dWW[j][i] = dWW[i][j];
1497  } }
1498 
1499  Di = toPars*dW;
1500 // Dij = (toPars*dWW)*(toPars.T());
1501  myTCL::mxmlrtS(toPars,dWW,Dij);
1502  for (int i=0;i<mNPars;i++) {
1503  for (int j=0;j<i;j++) {
1504  assert(fabs(dWW[i][j]-dWW[j][i])<1e-10);}}
1505  return Fcn;
1506 }
1507 //______________________________________________________________________________
1508 double HitPars_t::DERIV(int npt, const MyPull *pt,TVectorD &Di,TMatrixD &Dij,int maxTrk)
1509 {
1510 
1511  Di.ResizeTo(mNPars);
1512  Dij.ResizeTo(mNPars,mNPars);
1513  Di = 0.; Dij=0.; mNTrks=0;
1514  TVectorD Ti ,TiM (mNPars);
1515  TMatrixD Tij,TijM(mNPars,mNPars);
1516  int jl=0,trk=pt[0].trk;
1517  double xl=0;
1518  double sum=0,sumM=0;
1519  int NSave = (int)sqrt((npt/30.));
1520 
1521  for (int jr=1; 1; jr++) {
1522  double xl0 = xl; xl = pt[jr].xyz[0];
1523  if (jr<npt && pt[jr].trk==trk && xl0<xl) continue;
1524  int n = jr-jl;
1525  if (n>15) {
1526  double tsum = Deriv(n,pt+jl,Ti,Tij);
1527  sum += tsum;
1528  Di += Ti ;
1529  Dij += Tij;
1530  mNTrks++;
1531  if (mNTrks>=maxTrk) break;
1532  if (((mNTrks+1)%NSave)==0 ) {
1533  sumM += sum; sum= 0.;
1534  TiM += Di ; Di = 0.;
1535  TijM += Dij; Dij= 0.;
1536  }
1537  }
1538  if (jr>=npt) break;
1539  trk = pt[jr].trk;
1540  jl = jr;
1541  }
1542  sum +=sumM;
1543  Di +=TiM ;
1544  Dij +=TijM;
1545  return sum;
1546 }
1547 //______________________________________________________________________________
1548 void HitPars_t::myDers( double fake,double wy, double wz, double dens
1549  , double g[2][3],double G[2][3][3])
1550 {
1551 //A = 3.14/sqrt(wy*wz)
1552 //WY = wy*(1+2*fake*dens*A)
1553 //WZ = wz*(1+2*fake*dens*A)
1554  double u[3]={fake,wy,wz};
1555  double Cd[3],Cdd[3][3];
1556  double (*myCdd)[3] = Cdd;
1557  if (!G) myCdd=0;
1558  double C0 = myFake(fake,wy,wz,dens,Cd,myCdd);
1559 
1560  memset(g[0],0,2*3*sizeof(g[0][0]));
1561  for (int iyz=0;iyz<2;iyz++) {
1562  for (int ivar=0;ivar<3;ivar++) {
1563  g[iyz][ivar] = u[iyz+1]*Cd[ivar];
1564  if (ivar==iyz+1) g[iyz][ivar]+= C0; }}
1565 
1566  if (!G) return;
1567 
1568  memset(G[0][0],0,2*3*3*sizeof(G[0][0][0]));
1569  for (int iyz=0;iyz<2;iyz++) {
1570  for (int i=0;i<3;i++) {
1571  for (int j=0;j<3;j++) {
1572  G[iyz][i][j] = u[iyz+1]*Cdd[i][j];
1573  if (i==iyz+1) G[iyz][i][j]+=Cd[j];
1574  if (j==iyz+1) G[iyz][i][j]+=Cd[i];
1575  } } }
1576 
1577 }
1578 //______________________________________________________________________________
1579 double HitPars_t::myFake( double fake,double wy, double wz, double dens
1580  , double Cd[3],double Cdd[3][3])
1581 {
1582 //A = 3.14/sqrt(wy*wz)
1583 //WY = wy*(1+2*fake*dens*A)
1584 //WZ = wz*(1+2*fake*dens*A)
1585  double dens2 = 2*dens;
1586 
1587  double A = 3.14/sqrt(wy*wz);
1588  double Ay = -0.5*A/wy ;
1589  double Az = -0.5*A/wz ;
1590 
1591  double C0 = (1+fake*dens2*A);
1592  if (!Cd) return C0;
1593  Cd[0] = dens2*A;
1594  Cd[1] = fake*dens2*Ay;
1595  Cd[2] = fake*dens2*Az;
1596  if (!Cdd) return C0;
1597 
1598  double Ayy = 0.5*(-Ay + A/wy)/wy;
1599  double Ayz = 0.5*(-Az )/wy;
1600  double Azz = 0.5*(-Az + A/wz)/wz;
1601 
1602  Cdd[0][0] = 0;
1603  Cdd[0][1] = dens2*Ay;
1604  Cdd[0][2] = dens2*Az;
1605  Cdd[1][0] = Cdd[0][1];
1606  Cdd[1][1] = fake*dens2*Ayy;
1607  Cdd[1][2] = fake*dens2*Ayz;
1608  Cdd[2][0] = Cdd[0][2];
1609  Cdd[2][1] = Cdd[1][2];
1610  Cdd[2][2] = fake*dens2*Azz;
1611  return C0;
1612 }
1613 //______________________________________________________________________________
1614 void HitPars_t::Print(const HitPars_t *init) const
1615 {
1616  const char *TYZ[]={"y","z"};
1617  const char *bnd;
1618 
1619  printf("HitPars(Fake ) ");
1620  bnd = (Lim(0))? "*":" ";
1621  if (init) printf("\t ");
1622  printf("\tpout=%8.4g(+-%8.4g)%s\n",mDat[0],mDrr[0],bnd);
1623 
1624  for (int iDet=0;iDet<mNDets;iDet++) {
1625  for (int iYZ=0;iYZ<2;iYZ++) {
1626  int ln = Len(iDet,iYZ);
1627  int ipar0 = IPar(iDet,iYZ);
1628  assert(ipar0>=0);
1629  for (int ip=0;ip<ln;ip++) {
1630  bnd = (Lim(ipar0+ip))? "*":" ";
1631  double p = mPars[iDet][iYZ][ip];
1632  double e = mErrs[iDet][iYZ][ip];
1633  printf("HitPars(%d,%s,%d) ",iDet,TYZ[iYZ],ip);
1634  if (init) {
1635  printf("\tinit=%8.4g ",init->mPars[iDet][iYZ][ip]);}
1636  printf("\tpout=%8.4g(+-%8.4g)%s",p,e,bnd);
1637  printf("\n");
1638  } } }
1639 
1640 }
1641 //______________________________________________________________________________
1642 double HitPars_t::Diff(const HitPars_t &init) const
1643 {
1644  double pctMax=0;
1645  for (int i=0;i<mNPars;i++) {
1646  double dif = fabs((*this)[i]-init[i]);
1647  double err = init.Err(i);
1648  if (err<=0) err = 0.5*((*this)[i]+init[i])+1e-10;
1649  dif/= err;
1650  if (pctMax<dif) pctMax=dif;
1651  }
1652  return pctMax*100;
1653 }
1654 //______________________________________________________________________________
1655 int HitPars_t::Test(int npt, const MyPull *pt) const
1656 {
1657  enum {kNTrkTst=5};
1658  TVectorD Di ,DiT ;
1659  TMatrixD Dij,DijT;
1660  HitPars_t HP(*this);
1661  HP[0]=9;
1662  int npars = HP.NPars();
1663  double fcn = HP.DERIV(npt,pt,Di,Dij,kNTrkTst);
1664  for (int i=0;i<npars;i++) {
1665  double sav = HP[i];
1666  double delta = fabs(sav)*1e-3;
1667  if (delta<1e-6) delta=1e-6;
1668  HP[i]=sav+delta;
1669  double fcnT = HP.DERIV(npt,pt,DiT,DijT,kNTrkTst);
1670  HP[i]=sav;
1671  double ana = 0.5*(Di[i]+DiT[i]);
1672  double num = (fcnT-fcn)/delta;
1673  double dif = (num-ana)/(fabs(num+ana)+1e-5);
1674  if (fabs(dif)>0.001)
1675  printf("\nDer(%2d): ana = %g \tnum = %g \tdif = %g\n\n",i,ana,num,dif);
1676  for (int k=0;k<npars;k++) {
1677  ana = 0.5*(Dij[i][k]+DijT[i][k]);
1678  num = (DiT[k]-Di[k])/delta;
1679  dif = (num-ana)/(fabs(num+ana)+1e-5);
1680  if (fabs(dif)>0.001)
1681  printf("Der(%2d,%2d): ana = %g \tnum = %g \tdif = %g\n",i,k,ana,num,dif);
1682  }
1683 
1684  }
1685  return 0;
1686 }
1687 //______________________________________________________________________________
1688 double HitPars_t::Err(const double Pars[3],int nPars,const double A[3])
1689 {
1690 static const double tenMicrons = 1e-3;
1691 static const double min2Err = tenMicrons*tenMicrons;
1692 static const double max2Err = 1.;
1693  double err = TCL::vdot(A,Pars,nPars);
1694  if (nPars<=1) return err;
1695  if (err<min2Err) err=min2Err;
1696  if (err>max2Err) err=max2Err;
1697  return err;
1698 }
1699 //______________________________________________________________________________
1700 void HitPars_t::Limit()
1701 {
1702  for (int i=0;i<mNPars;i++) {
1703  if (mDat[i]<mMin[i]) mDat[i]=mMin[i];
1704  if (mDat[i]>mMax[i]) mDat[i]=mMax[i];
1705  }
1706 }
1707 #if 1
1708 //______________________________________________________________________________
1709 void HitPars_t::Prep(int npt, const MyPull *pt,TVectorD &Y,TVectorD &Z
1710  ,TVectorD &S,TVectorD &cos2Psi)
1711 {
1712 static int myDebug=0;
1713 
1714  Y.ResizeTo(npt);
1715  Z.ResizeTo(npt);
1716  S.ResizeTo(npt);
1717  if (myDebug) Show(npt,pt);
1718  int npt21 = npt*2+1;
1719  cos2Psi.ResizeTo(npt21);cos2Psi = 1.;
1720 
1721  THelixFitter helx;
1722  for (int ipt=0;ipt<npt;ipt++) {
1723  TVector3 point;
1724  point[0] = pt[ipt].xhg();
1725  point[1] = pt[ipt].yhg();
1726  point[2] = pt[ipt].zhg();
1727  helx.Add(point[0],point[1],point[2]);
1728  double emx[3]={pow(pt[ipt].hye,2),0,pow(pt[ipt].hye,2)};
1729  helx.AddErr(emx,pow(pt[ipt].hze,2));
1730  }
1731  helx.Fit();
1732  THelixTrack move(helx);
1733  double s = 0;
1734  for (int ipt=0;ipt<npt;ipt++) {
1735  TVector3 point;
1736  point[0] = pt[ipt].xhg();
1737  point[1] = pt[ipt].yhg();
1738  point[2] = pt[ipt].zhg();
1739  double ds = move.Path(point[0],point[1]);
1740  if (ipt) s+=ds;
1741  move.Move(ds);
1742  TVector3 pos(move.Pos());
1743  TVector3 dir(move.Dir());
1744  TVector3 nor(dir[1],-dir[0],0.); nor = nor.Unit();
1745  double dy = (point-pos)*nor;
1746  double dz = point[2]-pos[2];
1747  S[ipt]=s;
1748  Y[ipt]=dy;
1749  Z[ipt]=dz;
1750  cos2Psi[ipt+1]=pow(cos(pt[ipt].psi),2);
1751  }
1752 }
1753 #endif//1
1754 //______________________________________________________________________________
1755 void HitPars_t::Show(int npt, const MyPull *pt)
1756 {
1757 // lev=0 draw all nodes
1758 // lev=1 draw hit nodes
1759 // lev=2 draw hits only
1760 
1761 
1762 
1763 static TCanvas *myCanvas=0;
1764 static TGraph *graph[3][3] = {{0,0,0},{0,0,0},{0,0,0}};
1765 // P G P G
1766  float X[3][3][100],Y[3][3][100];
1767  int N[3][3];
1768 
1769  if (!myCanvas) {myCanvas = new TCanvas("C1","",600,800);
1770  myCanvas->Divide(1,3);}
1771  if(pt) {
1772  double curv = 0,xPrev,yPrev; int nCurv = 0;
1773  for (int i=0;i<9;i++) {delete graph[0][i];graph[0][i]=0;}
1774  for (int ig=0;ig<3;ig++) {
1775  int n=0;
1776  double s = 0;
1777  xPrev = 3e33;
1778  for (int ipt=0;ipt<npt;ipt++) {
1779  const MyPull *node = pt+ipt;
1780 // S calculation common based on node x,y for both hit and node
1781  double xNode = node->x_g();
1782  double yNode = node->y_g();
1783  if (xPrev<3e33) {
1784  double ds = sqrt(pow(xNode-xPrev,2)+pow(yNode-yPrev,2));
1785  double si = 0.5*ds*curv; if (si>0.99) si=0.99;
1786  if (si>0.01) ds = 2*asin(si)/curv;
1787  s += ds;
1788  }
1789  xPrev = xNode;
1790  yPrev = yNode;
1791 
1792  if (ig==0) { curv += node->curv; nCurv++; continue;}
1793  if (ig==1) {//draw nodes
1794  X[0][ig][n] = node->x_g();
1795  Y[0][ig][n] = node->y_g();
1796  Y[2][ig][n] = node->z_g();
1797  } else {//draw hits only
1798  X[0][ig][n] = node->xhg();
1799  Y[0][ig][n] = node->yhg();
1800  Y[2][ig][n] = node->zhg();
1801  }
1802 
1803  if (n) {
1804  float xh = X[0][ig][n]-X[0][ig][0];
1805  float yh = Y[0][ig][n]-Y[0][ig][0];
1806  float rh = xh*xh+yh*yh+1E-10;
1807  X[1][ig][n-1] = xh/rh;
1808  Y[1][ig][n-1] = yh/rh;
1809  }
1810  X[2][ig][n]=s;
1811  n++;
1812  }//end for nodes
1813  if (ig==0) { curv=fabs(curv)/nCurv; continue;}
1814  N[0][ig] = n;
1815  N[1][ig] = n-1;
1816  N[2][ig] = n;
1817  }//end for ig
1818 
1819  for (int ip=0;ip<3;ip++) {
1820  double xMin=999,xMax=-999,yMin=999,yMax=-999;
1821  for (int ig=1;ig<3;ig++) {
1822  for (int j=0;j<N[ip][ig];j++) {
1823  double x = X[ip][ig][j];
1824  if (xMin> x) xMin = x;
1825  if (xMax< x) xMax = x;
1826  double y = Y[ip][ig][j];
1827  if (yMin> y) yMin = y;
1828  if (yMax< y) yMax = y;
1829  }
1830  }
1831  X[ip][0][0] = xMin; Y[ip][0][0] = yMin;
1832  X[ip][0][1] = xMin; Y[ip][0][1] = yMax;
1833  X[ip][0][2] = xMax; Y[ip][0][2] = yMin;
1834  X[ip][0][3] = xMax; Y[ip][0][3] = yMax;
1835  N[ip][0] = 4;
1836  }
1837 static const char *opt[]={"AP","Same CP","Same *"};
1838  for (int ip=0;ip<3;ip++) {
1839  for (int ig =0;ig<3;ig++) {
1840  graph[ip][ig] = new TGraph(N[ip][ig] , X[ip][ig], Y[ip][ig]);
1841  if(ig==2) graph[ip][ig]->SetMarkerColor(kRed);
1842  myCanvas->cd(ip+1); graph[ip][ig]->Draw(opt[ig]);
1843  }//end for ig
1844  }//end ip
1845 
1846  }//end if
1847 
1848 
1849  if (!myCanvas) return;
1850  myCanvas->Modified();
1851  myCanvas->Update();
1852  while(!gSystem->ProcessEvents()){};
1853 }
1854 //______________________________________________________________________________
1855 int HitPars_t::Test()
1856 {
1857 return 0;
1858 }
1859 //______________________________________________________________________________
1860 //______________________________________________________________________________
1861 //______________________________________________________________________________
1862 //______________________________________________________________________________
1863 double myTCL::vmaxa(const double *a,int na)
1864 {
1865  double r=0;
1866  for (int i=0;i<na;i++){if (r < fabs(a[i])) r = fabs(a[i]);}
1867  return r;
1868 }
1869 //______________________________________________________________________________
1870 double myTCL::vmaxa(const TVectorD &a)
1871 {
1872  return vmaxa(a.GetMatrixArray(),a.GetNrows());
1873 }
1874 //______________________________________________________________________________
1875 void myTCL::vfill(double *a,double f,int na)
1876 {
1877  for (int i=0;i<na;i++) {a[i]=f;}
1878 }
1879 //_____________________________________________________________________________
1880 void myTCL::eigen2(const double err[3], double lam[2], double eig[2][2])
1881 {
1882 
1883  double spur = err[0]+err[2];
1884  double det = err[0]*err[2]-err[1]*err[1];
1885  double dis = spur*spur-4*det;
1886  if (dis<0) dis = 0;
1887  dis = sqrt(dis);
1888  lam[0] = 0.5*(spur+dis);
1889  lam[1] = 0.5*(spur-dis);
1890  eig[0][0] = 1; eig[0][1]=0;
1891  if (dis>1e-6*spur) {// eigenvalues are different
1892  if (fabs(err[0]-lam[0])>fabs(err[2]-lam[0])) {
1893  eig[0][1] = 1; eig[0][0]= -err[1]/(err[0]-lam[0]);
1894  } else {
1895  eig[0][0] = 1; eig[0][1]= -err[1]/(err[2]-lam[0]);
1896  }
1897  double tmp = sqrt(eig[0][0]*eig[0][0]+eig[0][1]*eig[0][1]);
1898  eig[0][0]/=tmp; eig[0][1]/=tmp;
1899  }
1900  eig[1][0]=-eig[0][1]; eig[1][1]= eig[0][0];
1901 }
1902 //______________________________________________________________________________
1903 /*
1904 * $Id: fiterr.C,v 1.8 2019/07/02 20:48:18 perev Exp $
1905 *
1906 * $Log: fiterr.C,v $
1907 * Revision 1.8 2019/07/02 20:48:18 perev
1908 * Increaese array of hits
1909 *
1910 * Revision 1.7 2010/01/27 21:38:06 perev
1911 * bugFix + remove not used code
1912 *
1913 * Revision 1.6 2008/01/20 00:43:02 perev
1914 * Mess with timestamps fixed
1915 *
1916 * Revision 1.5 2007/04/26 04:25:32 perev
1917 * Cleanup
1918 *
1919 * Revision 1.1.1.1 1996/04/01 15:02:13 mclareni
1920 * Mathlib gen
1921 */
1922 double myTCL::simpson(const double *F,double A,double B,int NP)
1923 {
1924  int N2=NP-1;
1925  assert(N2>0 && !(N2&1));
1926  double S1=F[N2-1];
1927  double S2=0;
1928 
1929  for (int N = 1;N<=N2-3;N+=2) {S1+=F[N];S2+=F[N+1];}
1930  S1=S1+S1+S2;
1931  double H=(F[0]+F[N2]+S1+S1)*(B-A)/(3*N2);
1932  return H;
1933 }
1934 //______________________________________________________________________________
1935 void myTCL::mxmlrt(const TMatrixD &A,const TMatrixD &B,TMatrixD &X)
1936 {
1937  int nRowA = A.GetNrows();
1938  int nColA = A.GetNcols();
1939  int nRowB = B.GetNrows();
1940  int nColB = B.GetNcols(); if(nColB){}
1941  assert(nColA ==nRowB);
1942  X.ResizeTo(nRowA,nRowA);
1943  TCL::mxmlrt(A.GetMatrixArray(),B.GetMatrixArray()
1944  ,X.GetMatrixArray(),nRowA,nColA);
1945 
1946 }
1947 //______________________________________________________________________________
1948 void myTCL::mxmlrtS(const double *A,const double *B,double *X,int nra,int nca)
1949 {
1950  TCL::vzero(X,nra*nra);
1951  for (int i=0,ii=0;i<nra;i++,ii+=nca) {
1952  for (int j=0,jj=0;j<nca;j++,jj+=nca) {
1953  if(!A[ii+j]) continue;
1954  for (int k=0,kk=0;k<=i;k++,kk+=nca) {
1955  double &Xik =X[i*nra+k];
1956  for (int l=0;l<nca;l++) {
1957  if(!A[kk+l]) continue;
1958  Xik +=A[ii+j]*A[kk+l]*B[jj+l];
1959  } } } }
1960  for (int i=0;i<nra;i++){
1961  for (int k=0;k<i ;k++){X[k*nra+i] = X[i*nra+k];}}
1962 }
1963 //______________________________________________________________________________
1964 void myTCL::mxmlrtS(const TMatrixD &A,const TMatrixD &B,TMatrixD &X)
1965 {
1966  int nRowA = A.GetNrows();
1967  int nColA = A.GetNcols();
1968  int nRowB = B.GetNrows();
1969  int nColB = B.GetNcols(); if(nColB){}
1970  assert(nColA ==nRowB);
1971  X.ResizeTo(nRowA,nRowA);
1972  myTCL::mxmlrtS(A.GetMatrixArray(),B.GetMatrixArray()
1973  ,X.GetMatrixArray(),nRowA,nColA);
1974 
1975 }
1976 //______________________________________________________________________________
1977 double myTCL::vasum(const double *a, int na)
1978 {
1979  double sum = 0;
1980  for (int i=0;i<na;i++) { sum += fabs(a[i]);}
1981  return sum;
1982 }
1983 //______________________________________________________________________________
1984 double myTCL::vasum(const TVectorD &a)
1985 {
1986  return vasum(a.GetMatrixArray(),a.GetNrows());
1987 }
1988 //______________________________________________________________________________
1989 int myTCL::SqProgSimple( TVectorD &x
1990  ,const TVectorD &g,const TMatrixD &G
1991  ,const TVectorD &Min
1992  ,const TVectorD &Max, int iAktp)
1993 {
1994 static int nCall=0; nCall++;
1995  enum {kINIT=0,kADDCONS,kFREECONS};
1996  int kase = kINIT;
1997  int nPars = g.GetNrows();
1998  TVectorD xx(x),gg(g),add(nPars);
1999  TArrayI Side(nPars);
2000  int nCons=0,addCons = -1,freCons=-1,freSide=0,addSide=0,con=0;
2001  double maxGra=3e33;
2002  int iAkt = iAktp;
2003  while(1946) {
2004 // Eliminate outdated constrains
2005  freCons=-1; freSide=0;
2006  if (nCons && kase==kFREECONS ) {
2007  double tryGra=kSMALL; freCons=-1;
2008  for (int ix=0;ix<nPars;ix++) {
2009  if(!Side[ix]) continue;
2010  double gra = gg[ix]*Side[ix];
2011  if (gra< tryGra) continue;
2012  if (gra>=maxGra) continue;
2013  freCons=ix; tryGra=gra;
2014  }
2015  if (freCons>=0) {
2016  maxGra = tryGra;
2017  freSide = Side[freCons];
2018  Side[freCons]=0;
2019  nCons--;
2020  }
2021  }
2022 
2023  if(kase==kFREECONS && freCons<0) { break;} //DONE
2024 
2025 // Make new matrix, etc...
2026  TMatrixD S(G);
2027  TVectorD B(gg);
2028  if (nCons) {
2029  for (int ix=0;ix<nPars;ix++) {
2030  if (Side[ix]==0) continue;
2031  for (int jx=0;jx<nPars;jx++) {S[ix][jx]=0; S[jx][ix]=0;}
2032  S[ix][ix]=1; B[ix]=0;
2033  } }
2034  if (iAkt==0 ) {
2035 
2036  double det=S.Determinant();
2037  if (fabs(det)<1e-100) return -99;
2038  S.Invert(0);
2039 // if (det<0) iAkt=1;
2040 // else add = (-1.)*(S*B);
2041  add = (-1.)*(S*B);
2042  double along = B*add;
2043  if (along>0) add*=(-1.);
2044  }
2045 
2046  if (iAkt==1 ) {
2047  double bb = (B*B);
2048  double bSb = (B*(S*B));
2049  double tau = -bb/(fabs(bSb)+3e-33);
2050  add = tau*B;
2051 
2052  }
2053  if(kase==kFREECONS && freSide) { //Free constrain case
2054  if (add[freCons]*freSide > -kSMALL) {
2055  Side[freCons]=freSide; nCons++; continue;}
2056  }
2057 
2058 // Do we need new constrain?
2059  double fak=1;
2060  addCons = -1; addSide = 0;
2061  con = 0;
2062  for (int ix=0;ix<nPars;ix++) {
2063  if (Side[ix]) {add[ix]=0; con = 100*con+ix+1;continue;}
2064  double xi = xx[ix]+fak*add[ix];
2065  if (xi < Min[ix]){fak = (Min[ix]-xx[ix])/add[ix]; addCons=ix;addSide=-1;}
2066  if (xi > Max[ix]){fak = (Max[ix]-xx[ix])/add[ix]; addCons=ix;addSide= 1;}
2067  assert(fak<=1. && fak>=0.);
2068  }
2069  add*=fak;
2070  xx+= add;
2071  gg += G*add;
2072  maxGra=3e33;
2073  kase = kFREECONS; if (!addSide) continue;
2074  kase = kADDCONS;
2075  xx[addCons] = (addSide<0)? Min[addCons]:Max[addCons];
2076 // Add new constrain;
2077  Side[addCons] = addSide ;nCons++;
2078 
2079  } //end of while(1)
2080  x = xx;
2081  return abs(con);
2082 }
2083 #endif
TTable::SavePrimitive
virtual void SavePrimitive(std::ostream &out, Option_t *option="")
Save a primitive as a C++ statement(s) on output stream &quot;out&quot;.
Definition: TTable.cxx:1810
HitAccr
Definition: fiterr.C:81
myTCL
Definition: fiterr.C:62
FitState_t
Definition: fiterr.C:682
MyPull
Definition: fiterr.C:90
TTable
Definition: TTable.h:48
HitPars_t
Definition: fiterr.C:125
Poli2
Definition: fiterr.C:194
Generated by   doxygen 1.8.5