#define __COMPILE__
#ifdef __COMPILE__
#include<iostream.h>
#include<TROOT.h>
#include<TFile.h>
#include<TF2.h>
#include<TVirtualFitter.h>
#include<TCanvas.h>
#include<TStopwatch.h>
#endif
//
//
//
void fcnk0(Int_t &npar, Double_t *gin, Double_t &f, Double_t *x, Int_t iflag)
{
    Double_t r01[3], r1, w1, f1, vp1, t1, hr1[3];
    Double_t r02[3], r2, w2, f2, vp2, t2, hr2[3];
    //
    //  t1,t2 free parameters of the helices
    //
    t1=x[0]; t2=x[1];
    Int_t k=1;
    for (Int_t i=0;i<3;i++) {++k;r01[i]=x[k];}
    ++k; r1=x[k];
    ++k; w1=x[k];
    ++k; f1=x[k];
    ++k; vp1=x[k];
    for (Int_t i=0;i<3;i++) {++k;r02[i]=x[k];}
    ++k; r2=x[k];
    ++k; w2=x[k];
    ++k; f2=x[k];
    ++k; vp2=x[k];
    hr1[0] = r01[0] + r1*TMath::Cos(w1*t1-f1);
    hr1[1] = r01[1] + r1*TMath::Sin(w1*t1-f1);
    hr1[2] = r01[2] + vp1 * t1; 
    hr2[0] = r02[0] + r2*TMath::Cos(w2*t2-f2);
    hr2[1] = r02[1] + r2*TMath::Sin(w2*t2-f2);
    hr2[2] = r02[2] + vp2 * t2; 
    Double_t distance = TMath::Sqrt( (hr1[0]-hr2[0])*(hr1[0]-hr2[0]) +
                                     (hr1[1]-hr2[1])*(hr1[1]-hr2[1]) +
				     (hr1[2]-hr2[2])*(hr1[2]-hr2[2]) );
    f=distance;
    return;
}

//
//  Do the job
//
void Helidebug() {
  UInt_t imin, idraw;
  Int_t par1, par2;
  Double_t r01[3], r02[3], r1, r2, w1, w2, f1, f2, vp1, vp2;
  Double_t t1min, t1max, t2min, t2max, lim;
  Double_t init1=-11., init2=10., errt1=1., errt2=1, t1ll=0, t1hl=0, t2ll=0, t2hl=0;

  //
  //  Read  First and Second Helix parameters
  //

  FILE *inf = fopen("helidist.inp","r");
  Int_t letto;
   letto=fscanf(inf,"%*10c %u ",&idraw);
   letto=fscanf(inf,"%*10c %lf ",&r01[0]);
   letto=fscanf(inf,"%*10c %lf ",&r01[1]);
   letto=fscanf(inf,"%*10c %lf ",&r01[2]);
   letto=fscanf(inf,"%*10c %lf ",&r1);
   letto=fscanf(inf,"%*10c %lf ",&w1);
   letto=fscanf(inf,"%*10c %lf ",&f1);
   letto=fscanf(inf,"%*10c %lf ",&vp1);
   
   letto=fscanf(inf,"%*10c %lf ",&r02[0]);
   letto=fscanf(inf,"%*10c %lf ",&r02[1]);
   letto=fscanf(inf,"%*10c %lf ",&r02[2]);
   letto=fscanf(inf,"%*10c %lf ",&r2);
   letto=fscanf(inf,"%*10c %lf ",&w2);
   letto=fscanf(inf,"%*10c %lf ",&f2);
   letto=fscanf(inf,"%*10c %lf ",&vp2);

   letto=fscanf(inf,"%*10c %u ",&imin);

   letto=fscanf(inf,"%*10c %lf ",&t1min);
   letto=fscanf(inf,"%*10c %lf ",&lim);

   letto=fscanf(inf,"%*10c %lf ",&init1);
   letto=fscanf(inf,"%*10c %lf ",&init2);
   letto=fscanf(inf,"%*10c %lf ",&errt1);
   letto=fscanf(inf,"%*10c %lf ",&errt2);
   letto=fscanf(inf,"%*10c %lf ",&t1ll);
   letto=fscanf(inf,"%*10c %lf ",&t1hl);
   letto=fscanf(inf,"%*10c %lf ",&t2ll);
   letto=fscanf(inf,"%*10c %lf ",&t2hl);
  fclose(inf);

  // phases are given as fraction of Pi
  f1*=TMath::Pi() ;
  f2*=TMath::Pi() ;

  //
  // Here you set the distance bi-dimensional plot limits
  //
  t1max= -t1min;
  t2min=t1min; t2max=t1max;

  //
  //   Minimization of the distance
  //
  
  Double_t arglist[100];
  TVirtualFitter *minuit = TVirtualFitter::Fitter(0,16);
  minuit->SetPrecision(1.0e-4); //default 1.e-6
  arglist[0]=2;
  minuit->ExecuteCommand("SET STRA", arglist, 1); //default 1
  printf("Starting timer\n");
  TStopwatch timer;
  timer.Start();
  minuit->SetFCN(fcnk0);
  //
  //  
  minuit->SetParameter(0, "t1",    init1, errt1, t1ll, t1hl);
  minuit->SetParameter(1, "t2",    init2, errt2, t2ll, t2hl);
  //
  minuit->SetParameter(2,  "r01x",  r01[0], 0.0001, 0, 0);
  minuit->SetParameter(3,  "r01y",  r01[1], 0.0001, 0, 0);
  minuit->SetParameter(4,  "r01z",  r01[2], 0.0001, 0, 0);
  minuit->SetParameter(5,  "r1",    r1,     0.0001, 0, 0);
  minuit->SetParameter(6,  "w1",    w1,     0.0001, 0, 0);
  minuit->SetParameter(7,  "f1",    f1,     0.0001, 0, 0);
  minuit->SetParameter(8,  "vp1",   vp1,    0.0001, 0, 0);
  minuit->SetParameter(9,  "r02x",  r02[0], 0.0001, 0, 0);
  minuit->SetParameter(10, "r02y",  r02[1], 0.0001, 0, 0);
  minuit->SetParameter(11, "r02z",  r02[2], 0.0001, 0, 0);
  minuit->SetParameter(12, "r2",    r2,     0.0001, 0, 0);
  minuit->SetParameter(13, "w2",    w2,     0.0001, 0, 0);
  minuit->SetParameter(14, "f2",    f2,     0.0001, 0, 0);
  minuit->SetParameter(15, "vp2",   vp2,    0.0001, 0, 0);
  
  minuit->FixParameter(2);
  minuit->FixParameter(3);
  minuit->FixParameter(4);
  minuit->FixParameter(5);
  minuit->FixParameter(6);
  minuit->FixParameter(7);
  minuit->FixParameter(8);
  minuit->FixParameter(9);
  minuit->FixParameter(10);
  minuit->FixParameter(11);
  minuit->FixParameter(12);
  minuit->FixParameter(13);
  minuit->FixParameter(14);
  minuit->FixParameter(15);
    
  arglist[0] = 0;
  minuit->ExecuteCommand("SET PRINT", arglist, 1);
  minuit->ExecuteCommand("MIGRAD", arglist, 0);
  minuit->ExecuteCommand("MINOS", arglist, 0);

  Double_t  finit1=0.0, ferrt1=0.0, ft1ll=0.0, ft1hl=0.0;
  Double_t  finit2=0.0, ferrt2=0.0, ft2ll=0.0, ft2hl=0.0; 
  //
  //
  // next line gives segmentation violation when compiled
  par1=minuit->GetParameter(0,"t1",finit1,ferrt1,ft1ll,ft1hl);
  par2=minuit->GetParameter(1,"t2",finit2,ferrt2,ft2ll,ft2hl);
  cout<<"t1 min= "<<finit1<<endl;
  cout<<"t2 min= "<<finit2<<endl;
  printf("Time at the end of job = %f seconds\n",timer.CpuTime());   
  return;
}