#include "/home/dbesson/src/gsl_multiroot_finder/Vtx4Hit.C"

void UVtxSpher4Hit(Bool_t LDebug0){
  Float_t xRx[MCH16], yRx[MCH16], zRx[MCH16], HitTime[MCH16];
  Float_t rZ, phiZ, theZ, dphiZ, dtheZ;
  Bool_t LUseAnt[MCH16]={false};
  Int_t nUseAntSpher=0;
  for(Int_t ich=0; ich<MCH16; ich++){
    if(LDebug)printf("Rx[%d]=%d HitTime=%g\n",ich,Rx[ich].LVtx,Rx[ich].HitTime);
    if(!Rx[ich].LVtx)continue;
    xRx[ich]=Rx[ich].x;
    yRx[ich]=Rx[ich].y;
    zRx[ich]=Rx[ich].z;
    LUseAnt[ich]=true;
    HitTime[ich]=Rx[ich].HitTime;
    //    HitTime[ich]=Rx[ich].HitTime;
    //23Nov17    LUseAnt[ich]=Rx[ich].iLVtx==1;
    if(LDebug)printf("Adding ch %d in UVtx4Hit\n",ich);
    if(xRx[ich]==0)printf("ERRZ: UVtx4Hit !!xRx[%d]=0\n",ich);
    if(yRx[ich]==0)printf("ERRZ: UVtx4Hit !!yRx[%d]=0\n",ich);
    if(zRx[ich]==0)printf("ERRZ: UVtx4Hit !!zRx[%d]=0\n",ich);
    if(HitTime[ich]==-NNN)printf("ERRZ!!!! HitTime[%d]=%g Rx.HitTime=%g Rx.LVtx=%d\n",ich,HitTime[ich],Rx[ich].HitTime,Rx[ich].LVtx);
    nUseAntSpher++;
  }
  printf("nUseAntSpher=%d\n",nUseAntSpher);
  Vtx4Hit(LDebug0,HitTime,xRx,yRx,zRx,LUseAnt,Fit4.r,Fit4.phi,Fit4.the,Fit4.dphi,Fit4.dthe,Fit4.ncomb,Fit4solnCh); //HitTime[0],&xRx[0],&yRx[0],&zRx[0],&LUseAnt[0]);
}

inline void FitMinuitXC(const char *rC, const char *phiC, const char *theC, Bool_t Lmnseek, Bool_t Lmigrad){
  TMinuit* eventrecoMinuit=new TMinuit(MVARMINUIT);
  eventrecoMinuit->SetPrintLevel(0); //-1 is minimum info 
  eventrecoMinuit->SetMaxIterations(MAXITERMINUIT);
  arglist[0]=200;
  arglist[1]=1;
  eventrecoMinuit->mnexcm("SET ERR",arglist,1,idum); // Parameter UP for both SIMPLEX and MIGRAD is set up here to 200
  //  eventrecoMinuit->mnrset(1);
  eventrecoMinuit->SetFCN(calc_chi_square_XC); //chi square routine - MUST RESET AFTER LOCAL CALL TO GAUS, E.G.!!!!!! CRASHES HERE FOR A1/ARAROOT 6 (6mar19)
  eventrecoMinuit->mnparm(0,rC,R0FITMINUIT,step[0],0,RMAXMINUIT,ierflg);
  eventrecoMinuit->mnparm(1,phiC,Rx[ichMinHitTime].phi,step[1],-PHIMAXMINUIT*2,PHIMAXMINUIT*2,ierflg);
  eventrecoMinuit->mnparm(2,theC,Rx[ichMinHitTime].the,step[2],0,THEMAXMINUIT*2,ierflg); 
  if(nFitVtx==3){
    eventrecoMinuit->mnparm(0,rC,1000.,step[0],0,RMAXMINUIT,ierflg);
    eventrecoMinuit->FixParameter(0);  
  }
  if(Lmnseek)eventrecoMinuit->mnseek(); //CPU-hog?
  eventrecoMinuit->mnsimp();
  if(Lmigrad)eventrecoMinuit->mnmigr(); 
  for(Int_t ivar=0; ivar<3; ivar++)eventrecoMinuit->GetParameter(ivar,FitPar[ivar],FitErr[ivar]);
  eventrecoMinuit->mnerrs(1,ddum,ddum1,FitErr[1],ddum2);
  eventrecoMinuit->mnerrs(2,ddum,ddum1,FitErr[2],ddum2);
  if(LMC)printf("eventrecoMinuit phi=%g the=%g dphi=%g dthe=%g\n",FitPar[1],FitPar[2],FitErr[1],FitErr[2]); 
  if(LMC)printf("Sim  r=%g    Simphi=%g            Simthe=%g\n",simVtx[0],simVtx[1],simVtx[2]);
  eventrecoMinuit->mnprin(3,fdum);
  printf("ierflg=%d / chisquare=%g convergence=%d\n",ierflg,ChiSqTot,convergence[0]); //ChiSqTot==0 =>failed fit!!!!
  printf("Fit4.r=%g Fit4.phi=%g+/-%g Fit4.the=%g+/-%g ncomb=%d\n",Fit4.r,Fit4.phi,Fit4.dphi,Fit4.the,Fit4.dthe,Fit4.ncomb);
  delete eventrecoMinuit;
}