// test Helix fit

doHelixA(int nEve=200) {
 gROOT->LoadMacro("$STAR/StRoot/StMuDSTMaker/COMMON/macros/loadSharedLibraries.C");
  loadSharedLibraries();
  cout << " loading done " << endl;


 assert( !gSystem->Load("ForwardHelix"));
 printf("ok\n");

  h1=new TH1F("chi2","chi2/DOF for helix fit",100,0.,5.);
  h2=new TH1F("rdPt","#Delta PT/PT == #Delta R/R",100,-2.,2.);
  h3=new TH1F("rdInvR","#Delta(1/R) / 1/R ;(1/Rgen - 1/Rreco)/ (1/Rgen) ",100,-2.,2.);
  h4=new TH2F("Dchi2","chi2/DOF vs. phi0; phi0 of prim track ; chi2/DOF",60,-3.14,3.14,100,0.,5.);


 //..... helix fit..............
  HelixUtil *fitter = new  HelixUtil;

  mRnd=new TRandom3;
  mRnd->SetSeed(0.0);
  printf(" First rnd=%f\n",mRnd->Uniform());

  const int mx=10;
  // all in cm
  double R=6600;
  double delPhi=0.020; // rad
  double sigArc=0.1; /*cm*/

  printf("Inp: R=%.2f lenArc/rad=%f, sigArc/cm=%.3f\n",R,delPhi,sigArc);

  double pi=2*acos(0);

  int iEve;
  //==============================================
  //  Events loop
  //==============================================
  for(iEve=0;iEve<nEve;iEve++) {

    double xA[mx],yA[mx],wei[mx];
    int nd=mx;
    
    // produce data points
    double W=1./sigArc/sigArc;
    int i;

    double phi0=(2*mRnd->Uniform()-1)*pi; // change direction of track over 12 sectors
    //  phi0=pi;

    //calc center of the circle
    double xc=R*cos(phi0);
    double yc=R*sin(phi0);
  
    printf("\nCenter          xc=%.2f yc=%.2f phi0/rad=%f\n",xc,yc,phi0);
    for(i=0;i<nd;i++) {
      
      double phi=-phi0+i*delPhi/nd +0.005; // measured from the center of the circle
      double xh=xc-R*cos(phi);
      double yh=yc+R*sin(phi);
      // add error along the padrow
      double er=mRnd->Gaus(0.,sigArc); // this is error in arc 
      double phiTpc=atan2(yh,xh);  // measured from the center of the TPC
      double ex=-er*sin(phiTpc);
      double ey=er*cos(phiTpc);
      //   ex=mRnd->Gaus(0.,sigArc);  ey=mRnd->Gaus(0.,sigArc); // 2D Gaussian

      xA[i]=xh+ex;
      yA[i]=yh+ey;
      wei[i]=W;
      //   if(i==0 || i==nd-1)
      // printf("IN: %d %.2f %.2f ex=%.2f ey=%.2f sq(W)=%.2f phiTpc=%.2f \n",i,xA[i],yA[i],ex,ey,sqrt(W),phiTpc);
    }
    
    double ddx=xA[0] - xA[nd-1];
    double ddy=yA[0] - yA[nd-1];
    double len=sqrt(ddx*ddx + ddy*ddy);

    fitter ->CircleFit(xA,yA,wei,nd);
    double rR=fitter->mFitRadius;
    double rxc=fitter->mXCenter;
    double ryc=fitter->mYCenter;
    double chi2=fitter-> mChi2Rad;
    double chi2dof=chi2/(1.*nd-3);
    int nIter=fitter-> mIterRad;
        
    printf("Fit: R=%.2f  xc=%.2f yc=%.2f chi2dof=%g nIter=%d len/cm=%.1f\n",rR,rxc,ryc,chi2dof,nIter,len);
    printf("DEL: R=%.2f  xc=%.2f yc=%.2f\n",rR-R,rxc-xc,ryc-yc);
    //  printf("D/R: R=%.2f  xc=%.2f yc=%.2f\n",(rR-R)/R,(rxc-xc)/R,(ryc-yc)/R);
    
    h1->Fill(chi2dof);
    h2->Fill((rR-R)/R);
    h3->Fill((1./rR-1./R)*R);
    h4->Fill(phi0,chi2dof);

    //======================= End of events ==================    
  }
}