--- DarthVader/TrackerLevel2/src/ExtTrack.cpp	2014/02/27 11:24:43	1.1
+++ DarthVader/TrackerLevel2/src/ExtTrack.cpp	2014/08/07 16:04:14	1.5
@@ -6,200 +6,591 @@
 #include <iostream>
 #include <math.h>
 using namespace std;
-//......................................
-// F77 routines
-//......................................
-extern "C" {    
-//     void dotrack_(int*, double*, double*, double*, double*, int*);
-//     void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
-//     void dotrack3_(int*, double*, double*, double*, double*,double*, double*, double*,double*,int*);
-//     void mini2_(int*,int*,int*); 
-//     void guess_();
-//     void gufld_(float*, float*);
-//     float risxeta2_(float *);
-//     float risxeta3_(float *);
-//     float risxeta4_(float *);
-//     float risyeta2_(float *);
-}
 
+ExtTrack::ExtTrack(int dim){
+
+  //cout << " ExtTrack::ExtTrack("<<dim<<")"<<endl;
+  xgood   = NULL;
+  ygood   = NULL;
+  multmaxx   = NULL;
+  multmaxy   = NULL;
+  xm      = NULL;
+  ym      = NULL;
+  zm      = NULL;
+  xma     = NULL;
+  yma     = NULL;
+  zma     = NULL;
+  xmb     = NULL;
+  ymb     = NULL;
+  zmb     = NULL;
+  resx    = NULL;
+  resy    = NULL;
+  xv      = NULL;
+  yv      = NULL;
+  zv      = NULL;
+  axv     = NULL;
+  ayv     = NULL;
+  dedx_x  = NULL;
+  dedx_y  = NULL;
+
+  // xGF = new float[TrkParams::nGF];
+  // yGF = new float[TrkParams::nGF];
+  // init vectors  
+  SetDimension(dim);//allocate vectors
+  // init fit variables
+  Reset();
+  // set fit parameters (default)
+  SetZ0(23.5);        
+  SetMiniDefault();  
+};
 //--------------------------------------
 //
 //
 //--------------------------------------
-ExtTrack::ExtTrack(){
-    chi2  = 0;
-    nstep = 0;
-    for(int it1=0;it1<5;it1++){
-	al[it1] = 0;
-	for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
-    };
-    nplanes=0;
-
-    xgood   = 0;
-    ygood   = 0;
-    xm      = 0;
-    ym      = 0;
-    zm      = 0;
-    resx    = 0;
-    resy    = 0;
-    xv      = 0;
-    yv      = 0;
-    zv      = 0;
-    axv     = 0;
-    ayv     = 0;
-    dedx_x  = 0;
-    dedx_y  = 0;
-
-    TrkParams::SetMiniDefault();
+ExtTrack::ExtTrack( const ExtTrack& t ){
+  //cout << " ExtTrack::ExtTrack( const ExtTrack& t ) "<<this<<endl;
 
+  // --------------------
+  // initialize pointers
+  // --------------------
+
+  xgood   = NULL;
+  ygood   = NULL;
+  multmaxx   = NULL;
+  multmaxy   = NULL;
+  xm      = NULL;
+  ym      = NULL;
+  zm      = NULL;
+  xma     = NULL;
+  yma     = NULL;
+  zma     = NULL;
+  xmb     = NULL;
+  ymb     = NULL;
+  zmb     = NULL;
+  resx    = NULL;
+  resy    = NULL;
+  xv      = NULL;
+  yv      = NULL;
+  zv      = NULL;
+  axv     = NULL;
+  ayv     = NULL;
+  dedx_x  = NULL;
+  dedx_y  = NULL;
+
+  // xGF = new float[TrkParams::nGF];
+  // yGF = new float[TrkParams::nGF];
+  
+  SetDimension(t.nplanes);//allocate vectors
+
+  // --------------------
+  // copy values
+  // --------------------
+  zini  = t.zini;
+  chi2  = t.chi2;
+  nstep = t.nstep;
+  for(int it1=0;it1<5;it1++){
+    al[it1] = t.al[it1];
+    for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
+  };
+
+
+  for(int ip=0;ip<nplanes;ip++){
+    xgood[ip]  = t.xgood[ip];
+    ygood[ip]  = t.ygood[ip];
+    multmaxx[ip]  = t.multmaxx[ip];
+    multmaxy[ip]  = t.multmaxy[ip];
+    xm[ip]     = t.xm[ip];
+    ym[ip]     = t.ym[ip];
+    zm[ip]     = t.zm[ip];
+    xma[ip]     = t.xma[ip];
+    yma[ip]     = t.yma[ip];
+    zma[ip]     = t.zma[ip];
+    xmb[ip]     = t.xmb[ip];
+    ymb[ip]     = t.ymb[ip];
+    zmb[ip]     = t.zmb[ip];
+    resx[ip]   = t.resx[ip];
+    resy[ip]   = t.resy[ip];
+    xv[ip]     = t.xv[ip];
+    yv[ip]     = t.yv[ip];
+    zv[ip]     = t.zv[ip];
+    axv[ip]    = t.axv[ip];
+    ayv[ip]    = t.ayv[ip];
+    dedx_x[ip] = t.dedx_x[ip];
+    dedx_y[ip] = t.dedx_y[ip];
+  };
+  
+  int ngf = TrkParams::nGF;
+  for(int i=0; i<ngf; i++){
+    xGF[i] = t.xGF[i];
+    yGF[i] = t.yGF[i];
+  }
+  
 };
 //--------------------------------------
 //
 //
 //--------------------------------------
-ExtTrack::ExtTrack(Int_t dim){
-    chi2  = 0;
-    nstep = 0;
-    for(int it1=0;it1<5;it1++){
-	al[it1] = 0;
-	for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
-    };
-    
-    SetDimension(dim);
+void ExtTrack::Copy( ExtTrack& t ){
+
+  //  cout << " ExtTrack::Copy( ExtTrack& t ) "<<endl;
 
-    TrkParams::SetMiniDefault();
+  t.Reset();//reset variables
+  t.Clear();//deallocate vectors
+  t.SetDimension(nplanes);//allocate vectors
+
+  t.zini  = zini;
+  t.chi2  = chi2;
+  t.nstep = nstep;
+  for(int it1=0;it1<5;it1++){
+    t.al[it1] = al[it1];
+    for(int it2=0;it2<5;it2++)t.coval[it1][it2] = coval[it1][it2];
+  };    
+  //    cout << "ExtTrack::Copy()........nplanes "<<nplanes<<endl;
+  for(int ip=0;ip<nplanes;ip++){
+    t.xgood[ip]  = xgood[ip];
+    t.ygood[ip]  = ygood[ip];
+    t.multmaxx[ip]  = multmaxx[ip];
+    t.multmaxy[ip]  = multmaxy[ip];
+    t.xm[ip]     = xm[ip];
+    t.ym[ip]     = ym[ip];
+    t.zm[ip]     = zm[ip];
+    t.xma[ip]     = xma[ip];
+    t.yma[ip]     = yma[ip];
+    t.zma[ip]     = zma[ip];
+    t.xmb[ip]     = xmb[ip];
+    t.ymb[ip]     = ymb[ip];
+    t.zmb[ip]     = zmb[ip];
+    t.resx[ip]   = resx[ip];
+    t.resy[ip]   = resy[ip];
+    t.xv[ip]     = xv[ip];
+    t.yv[ip]     = yv[ip];
+    t.zv[ip]     = zv[ip];
+    t.axv[ip]    = axv[ip];
+    t.ayv[ip]    = ayv[ip];
+    t.dedx_x[ip] = dedx_x[ip];
+    t.dedx_y[ip] = dedx_y[ip];
+  };
+  int ngf = TrkParams::nGF;
+  for(int i=0; i<ngf; i++){
+    t.xGF[i] = xGF[i];
+    t.yGF[i] = yGF[i];
+  }
+  
 };
 //--------------------------------------
 //
 //
 //--------------------------------------
-ExtTrack::ExtTrack(const ExtTrack& t){
-    chi2  = t.chi2;
-    nstep = t.nstep;
-    for(int it1=0;it1<5;it1++){
-	al[it1] = t.al[it1];
-	for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
-    };
-    SetDimension(t.nplanes);
-    for(int ip=0;ip<nplanes;ip++){
-	xgood[ip]  = t.xgood[ip];
-	ygood[ip]  = t.ygood[ip];
-	xm[ip]     = t.xm[ip];
-	ym[ip]     = t.ym[ip];
-	zm[ip]     = t.zm[ip];
-	resx[ip]   = t.resx[ip];
-	resy[ip]   = t.resy[ip];
-	xv[ip]     = t.xv[ip];
-	yv[ip]     = t.yv[ip];
-	zv[ip]     = t.zv[ip];
-	axv[ip]    = t.axv[ip];
-	ayv[ip]    = t.ayv[ip];
-	dedx_x[ip] = t.dedx_x[ip];
-	dedx_y[ip] = t.dedx_y[ip];
-    };
-
-    TrkParams::SetMiniDefault();
+void ExtTrack::Set( TrkTrack& t , int index){
+
+  //  cout << "Set("<< &t <<", "<<index<<")"<<endl;
+  if(index < 0 ){ //NON FUNZIONA
+    cout << "Clear() "<<endl; 
+    Reset();//reset variables
+    Clear();//deallocate vectors
+    cout << "SetDimension("<<t.GetNhit()<<")"<<endl; 
+    SetDimension(t.GetNhit()); //allocate vectors
+//   }else{
+//       SetDimension(6+index);
+  }
+  // default fit settings (the same as TrkTrack)
+  SetMiniDefault();
+  SetZ0(23.5); 
+  //
+  int ipp;
+  if(index<0) ipp = 0;
+  else        ipp = index;
+  //
+//  cout << "Assign... nplanes = "<<nplanes<<endl;
+  for(int ip=0;ip<6;ip++){
+    
+    if(ipp>=nplanes){
+      cout << "void ExtTrack::Set( TrkTrack& , "<<index<<") -- ipp="<<ipp<<" exceed vector dimention "<<nplanes<<endl;
+      continue;
+    }
+
+    if( t.XGood(ip) || t.YGood(ip) || !(index<0) ){
+
+//      cout << "ipp "<<ipp<<" ip "<<ip<<endl;
+
+
+      xgood[ipp]  = t.XGood(ip); //NB! 
+      ygood[ipp]  = t.YGood(ip); //NB! 
+      multmaxx[ipp]  = t.multmaxx[ip]; 
+      multmaxy[ipp]  = t.multmaxy[ip]; 
+      xm[ipp]     = t.xm[ip];
+      ym[ipp]     = t.ym[ip];
+      zm[ipp]     = t.zm[ip];
+
+      float dummy = -100.;
+      xma[ipp]     = dummy;
+      yma[ipp]     = dummy;
+      zma[ipp]     = zm[ipp];
+      xmb[ipp]     = dummy;
+      ymb[ipp]     = dummy;
+      zmb[ipp]     = zm[ipp];
+
+      if( t.XGood(ip)*t.YGood(ip) ){ //double hit
+	xma[ipp]     = t.xm[ip];
+	yma[ipp]     = t.ym[ip];
+	zma[ipp]     = t.zm[ip];
+	xmb[ipp]     = t.xm[ip];
+	ymb[ipp]     = t.ym[ip];
+	zmb[ipp]     = t.zm[ip];
+      } else if (t.XGood(ip) || t.YGood(ip) ){                         //single hit
+	TrkParams::Load(5); 
+	if( !TrkParams::IsLoaded(5) ){
+	  cout << "void ExtTrack::FillMiniStruct(TrkTrack&) ---ERROR--- trk align.param. not loaded "<<endl;
+	  return;
+	}
+	double segment = 7.;// 2.;//cm //Elena 10th
+	int is = (int)t.GetSensor(ip); if(ip==5)is=1-is;
+	int ippp = 5-ip;
+	int il = (int)t.GetLadder(ip);	    
+	// NB: i parametri di allineamento hanno una notazione particolare!!!
+	// sensor = 0 (hybrid side), 1 
+	// ladder = 0-2 (increasing x)
+	// plane  = 0-5 (from bottom to top!!!) 
+	double omega   = 0.; 
+	double beta   = 0.; 
+	double gamma   = 0.; 
+	if( 
+	   (is < 0 || is > 1 || ip < 0 || ippp > 5 || il < 0 || il > 2) &&
+	   true){
+	  cout << " void ExtTrack::FillMiniStruct(TrkTrack&) --- WARNING ---trk sensor not defined, cannot read alignment parameters "<<endl;
+	  cout << " is ip il = "<<is<<" "<<ippp<<" "<<il<<endl;
+	}else{
+	  omega   = alignparameters_.omega[is][il][ippp];
+	  beta    = alignparameters_.beta[is][il][ippp];
+	  gamma   = alignparameters_.gamma[is][il][ippp];
+	}	  
+	if(       t.XGood(ip) && !t.YGood(ip) ){
+	  xma[ipp] = t.xm[ip] - omega * segment;
+	  yma[ipp] = t.ym[ip] + segment;
+	  zma[ipp] = t.zm[ip] + beta * segment;//not used yet
+	  xmb[ipp] = t.xm[ip] + omega * segment;
+	  ymb[ipp] = t.ym[ip] - segment;
+	  zmb[ipp] = t.zm[ip] - beta * segment;//not used yet
+	}else if( !t.XGood(ip) && t.YGood(ip) ){
+	  xma[ipp] = t.xm[ip] + segment;
+	  yma[ipp] = t.ym[ip] + omega * segment;
+	  zma[ipp] = t.zm[ip] - gamma * segment;//not used yet
+	  xmb[ipp] = t.xm[ip] - segment;
+	  ymb[ipp] = t.ym[ip] - omega * segment;
+	  zmb[ipp] = t.zm[ip] + gamma * segment;//not used yet
+	}
+      }
+	
+      // tailx[ipp]  = t.tailx[ip];
+      // taily[ipp]  = t.taily[ip];
+      resx[ipp]   = t.resx[ip];
+      resy[ipp]   = t.resy[ip];
+      xv[ipp]     = t.xv[ip];
+      yv[ipp]     = t.yv[ip];
+      zv[ipp]     = t.zv[ip];
+      axv[ipp]    = t.axv[ip];
+      ayv[ipp]    = t.ayv[ip];
+      dedx_x[ipp] = t.dedx_x[ip];
+      dedx_y[ipp] = t.dedx_y[ip];
+
+      ipp++;
+
+    }
+  };//endl loop on planes
+
+  chi2  = t.chi2;
+  nstep = t.nstep;
+  for(int it1=0;it1<5;it1++){
+    al[it1] = t.al[it1];
+    for(int it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
+  }; 
+ 
 };
 //--------------------------------------
 //
 //
 //--------------------------------------
-void ExtTrack::SetDimension(Int_t dim){
-
-    nplanes = dim;
-
-    xgood   = new Int_t[nplanes];
-    ygood   = new Int_t[nplanes];
-    xm      = new Float_t[nplanes];
-    ym      = new Float_t[nplanes];
-    zm      = new Float_t[nplanes];
-    resx    = new Float_t[nplanes];
-    resy    = new Float_t[nplanes];
-    xv      = new Float_t[nplanes];
-    yv      = new Float_t[nplanes];
-    zv      = new Float_t[nplanes];
-    axv     = new Float_t[nplanes];
-    ayv     = new Float_t[nplanes];
-    dedx_x  = new Float_t[nplanes];
-    dedx_y  = new Float_t[nplanes];
-    
-    SetMeasure(0,0,0);
-    SetResolution(0,0);
-    SetGood(0,0);
-    FitReset();
+void ExtTrack::SetMiniDefault(){
+  //  cout << " ExtTrack::SetMiniDefault()"<<endl;
+  exttrack_.trackmode = TrkParams::init__mini_trackmode;
+  exttrack_.fact      = TrkParams::init__mini_fact;   
+  exttrack_.istepmin  = TrkParams::init__mini_istepmin; 
+  TrkParams::SetDeltaB();
+  TrkParams::SetDLT();
 }
 //--------------------------------------
 //
 //
 //--------------------------------------
-void ExtTrack::Clear(Option_t* option){
-
-    FitReset();
-    if(nplanes){
-	delete []	xgood;
-	delete []	ygood;
-	delete []	xm;
-	delete []	ym;
-	delete []	zm;
-	delete []	resx;
-	delete []	resy;
-	delete []	xv;
-	delete []	yv;
-	delete []	zv;
-	delete []	axv;
-	delete []	ayv;
-	delete []	dedx_x;
-	delete []	dedx_y;
-    }
+void ExtTrack::SetDimension(int dim){
 
+  //cout << " ExtTrack::SetDimension("<<dim<<")"<<endl;;
+  if(dim<=0){
     nplanes = 0;
-    xgood   = 0;
-    ygood   = 0;
-    xm      = 0;
-    ym      = 0;
-    zm      = 0;
-    resx    = 0;
-    resy    = 0;
-    xv      = 0;
-    yv      = 0;
-    zv      = 0;
-    axv     = 0;
-    ayv     = 0;
-    dedx_x  = 0;
-    dedx_y  = 0;
+    if(dim<0)cout << "void ExtTrack::SetDimension("<<dim<<") --- invalid dimension "<<endl;
+    return;
+  }
+
+  //  Clear();
+
+  nplanes = dim;
+  
+  xgood   = new int[nplanes];
+  ygood   = new int[nplanes];
+  multmaxx   = new int[nplanes];
+  multmaxy   = new int[nplanes];
+  xm      = new float[nplanes];
+  ym      = new float[nplanes];
+  zm      = new float[nplanes];
+  xma     = new float[nplanes];
+  yma     = new float[nplanes];
+  zma     = new float[nplanes];
+  xmb     = new float[nplanes];
+  ymb     = new float[nplanes];
+  zmb     = new float[nplanes];
+  resx    = new float[nplanes];
+  resy    = new float[nplanes];
+  xv      = new float[nplanes];
+  yv      = new float[nplanes];
+  zv      = new float[nplanes];
+  axv     = new float[nplanes];
+  ayv     = new float[nplanes];
+  dedx_x  = new float[nplanes];
+  dedx_y  = new float[nplanes];
+  
+
+}
+//--------------------------------------
+//
+//
+//--------------------------------------
+bool ExtTrack::SetZ(int ip,float zmeas){
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::SetZ("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return false;
+  }
+  zm[ip] = zmeas;
+  zma[ip] = zmeas;
+  zmb[ip] = zmeas;
+  return true;
+
+};
+
+bool ExtTrack::ResetXY(int ip){
+
+  // cout << " ExtTrack::ResetXY("<<ip<<") "<<endl;
+
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::ResetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return false;
+  } 
+  xm[ip] = -100.;
+  ym[ip] = -100.;
+  xma[ip] = -100.;
+  yma[ip] = -100.;
+  xmb[ip] = -100.;
+  ymb[ip] = -100.;
+  resx[ip] = 1000.;
+  resy[ip] = 1000.;
+  dedx_x[ip] = 0.;
+  dedx_y[ip] = 0.;
+  xgood[ip] = 0;
+  ygood[ip] = 0;
+  multmaxx[ip] = 0;
+  multmaxy[ip] = 0;
+
+  return true;
+
+};
+
+bool ExtTrack::SetXY(int ip,float xmeas, float ymeas, float rx, float ry){
 
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::SetXY("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return false;
+  }
+  
+  xm[ip] = xmeas;
+  ym[ip] = ymeas;
+  xma[ip] = xmeas;
+  yma[ip] = ymeas;
+  xmb[ip] = xmeas;
+  ymb[ip] = ymeas;
+  resx[ip] = rx;
+  resy[ip] = ry;
+  // dedx_x[ip] = dedxx;
+  // dedx_y[ip] = dedxy;
+  xgood[ip] = 1;
+  ygood[ip] = 1;
+
+  return true;
+
+}; 
+bool ExtTrack::SetX(int ip,float xa, float xb, float ya, float yb,float res ){
+
+  // cout << " ExtTrack::SetX("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
+
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return false;
+  }
+  xm[ip] = -100.;
+  ym[ip] = -100.;
+  xma[ip] = xa;
+  yma[ip] = ya;
+  xmb[ip] = xb;
+  ymb[ip] = yb;
+  resx[ip] = res;
+  resy[ip] = 1000.;
+  xgood[ip] = 1;
+  ygood[ip] = 0;
+  // dedx_x[ip] = dedx;
+  // dedx_y[ip] = 0.;
+  return true;
+}; 
+bool ExtTrack::SetY(int ip,float xa, float xb, float ya, float yb,float res ){
+
+  // cout << " ExtTrack::SetY("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
+
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::SetX("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return false;
+  }
+  xm[ip] = -100.;
+  ym[ip] = -100.;
+  xma[ip] = xa;
+  yma[ip] = ya;
+  xmb[ip] = xb;
+  ymb[ip] = yb;
+  resx[ip] = 1000.;
+  resy[ip] = res;
+  xgood[ip] = 0;
+  ygood[ip] = 1;
+  // dedx_x[ip] = 0.;
+  // dedx_y[ip] = dedx;
+  return true;
+}; 
+
+bool ExtTrack::SetXGood(int ip,int icl_piu_uno, int il, int is ){
+
+  // cout << " ExtTrack::SetXGood("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
+
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::SetXGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return false;
+  }
+  xgood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno; 
+  return true;
+}; 
+bool ExtTrack::SetYGood(int ip,int icl_piu_uno, int il, int is ){ 
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::SetYGood("<<ip<<",...)  -- nplanes"<<nplanes<<endl;
+    cout << " ExtTrack::SetYGood("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return false;
+  }
+  ygood[ip]=(il+1)*100000000+(is+1)*10000000+icl_piu_uno; 
+  return true;
 };
+
 //--------------------------------------
 //
 //
 //--------------------------------------
-void ExtTrack::Delete(){
-    Clear();
+void ExtTrack::Clear(Option_t* option){
+
+  //  cout << " ExtTrack::Clear("<<option<<")  "<<this<<endl;
+
+
+  //  cout << " xgood "<<xgood<<endl;
+  
+  //  if(nplanes>0){
+  if(xgood) delete []	xgood;
+  if(ygood) delete []	ygood;
+  if(multmaxx) delete [] multmaxx;
+  if(multmaxy) delete [] multmaxy;
+  if(xm)    delete []	xm;
+  if(ym)    delete []	ym;
+  if(zm)    delete []	zm;
+  if(xma)   delete []	xma;
+  if(yma)   delete []	yma;
+  if(zma)   delete []	zma;
+  if(xmb)   delete []	xmb;
+  if(ymb)   delete []	ymb;
+  if(zmb)   delete []	zmb;
+  if(resx)  delete []	resx;
+  if(resy)  delete []	resy;
+  if(xv)    delete []	xv;
+  if(yv)    delete []	yv;
+  if(zv)    delete []	zv;
+  if(axv)   delete []	axv;
+  if(ayv)   delete []	ayv;
+  if(dedx_x)delete []	dedx_x;
+  if(dedx_y)delete []	dedx_y;
+  //}
+
+  xgood   = NULL;
+  ygood   = NULL;
+  multmaxx   = NULL;
+  multmaxy   = NULL;
+  xm      = NULL;
+  ym      = NULL;
+  zm      = NULL;
+  xma     = NULL;
+  yma     = NULL;
+  zma     = NULL;
+  xmb     = NULL;
+  ymb     = NULL;
+  zmb     = NULL;
+  resx    = NULL;
+  resy    = NULL;
+  xv      = NULL;
+  yv      = NULL;
+  zv      = NULL;
+  axv     = NULL;
+  ayv     = NULL;
+  dedx_x  = NULL;
+  dedx_y  = NULL;
+
+  nplanes = 0;
+
+  //  Reset();
+
 };
 
-/**
- * Set the  position measurements 
- */
-void ExtTrack::SetMeasure(Double_t *xmeas, Double_t *ymeas, Double_t *zmeas){
-    for(Int_t i=0; i<nplanes; i++) xm[i]= (xmeas? *xmeas++ : 0.);
-    for(Int_t i=0; i<nplanes; i++) ym[i]= (ymeas? *ymeas++ : 0.);
-    for(Int_t i=0; i<nplanes; i++) zm[i]= (ymeas? *zmeas++ : 0.);
-}
-/**
- * Set the TrkTrack position resolution 
- */
-void ExtTrack::SetResolution(Double_t *rx, Double_t *ry){
-    for(Int_t i=0; i<nplanes; i++) resx[i]= (rx? *rx++ : 0.);
-    for(Int_t i=0; i<nplanes; i++) resy[i]= (ry? *ry++ : 0.);
-}
-/**
- * Set the TrkTrack good measurement 
- */
-void ExtTrack::SetGood(Int_t *xg, Int_t *yg){
+void ExtTrack::Delete(){ 
+  Clear();
+  // delete [] xGF; 
+  // delete [] yGF; 
+}        
 
-    for(Int_t i=0; i<nplanes; i++) xgood[i]=(xg? *xg++ : 0 );
-    for(Int_t i=0; i<nplanes; i++) ygood[i]=(yg? *yg++ : 0 );
-}
+//--------------------------------------
+//
+//
+//--------------------------------------
+
+// /**
+//  * Set the  position measurements 
+//  */
+// void ExtTrack::SetMeasure(Double_t *xmeas, Double_t *ymeas, Double_t *zmeas){
+//     for(int i=0; i<nplanes; i++) xm[i]= (xmeas? *xmeas++ : 0.);
+//     for(int i=0; i<nplanes; i++) ym[i]= (ymeas? *ymeas++ : 0.);
+//     for(int i=0; i<nplanes; i++) zm[i]= (ymeas? *zmeas++ : 0.);
+// }
+// /**
+//  * Set the TrkTrack position resolution 
+//  */
+// void ExtTrack::SetResolution(Double_t *rx, Double_t *ry){
+//     for(int i=0; i<nplanes; i++) resx[i]= (rx? *rx++ : 0.);
+//     for(int i=0; i<nplanes; i++) resy[i]= (ry? *ry++ : 0.);
+// }
+// /**
+//  * Set the TrkTrack good measurement 
+//  */
+// void ExtTrack::SetGood(int *xg, int *yg){
+
+//     for(int i=0; i<nplanes; i++) xgood[i]=(xg? *xg++ : 0 );
+//     for(int i=0; i<nplanes; i++) ygood[i]=(yg? *yg++ : 0 );
+// }
 
 //--------------------------------------
 //
@@ -207,6 +598,7 @@
 //--------------------------------------
 void ExtTrack::Dump(){
     cout << endl << "========== Track " ;
+    cout << endl << "zini     : "<< zini;
     cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " ";
     cout << endl << "chi^2    : "<< chi2;
     cout << endl << "n.step   : "<< nstep;
@@ -215,6 +607,12 @@
     cout << endl << "xm       : "; for(int i=0; i<nplanes; i++)cout << xm[i] << " ";
     cout << endl << "ym       : "; for(int i=0; i<nplanes; i++)cout << ym[i] << " ";
     cout << endl << "zm       : "; for(int i=0; i<nplanes; i++)cout << zm[i] << " ";
+    cout << endl << "xma       : "; for(int i=0; i<nplanes; i++)cout << xma[i] << " ";
+    cout << endl << "yma       : "; for(int i=0; i<nplanes; i++)cout << yma[i] << " ";
+    cout << endl << "zma       : "; for(int i=0; i<nplanes; i++)cout << zma[i] << " ";
+    cout << endl << "xmb       : "; for(int i=0; i<nplanes; i++)cout << xmb[i] << " ";
+    cout << endl << "ymb       : "; for(int i=0; i<nplanes; i++)cout << ymb[i] << " ";
+    cout << endl << "zmb       : "; for(int i=0; i<nplanes; i++)cout << zmb[i] << " ";
     cout << endl << "xv       : "; for(int i=0; i<nplanes; i++)cout << xv[i] << " ";
     cout << endl << "yv       : "; for(int i=0; i<nplanes; i++)cout << yv[i] << " ";
     cout << endl << "zv       : "; for(int i=0; i<nplanes; i++)cout << zv[i] << " ";
@@ -225,90 +623,366 @@
     cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[2][i]<<" ";
     cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[3][i]<<" ";
     cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[4][i]<<" ";
+    cout << endl << "MDR      : "<< (coval[4][4]>0. ?  1./sqrt(coval[4][4]) : 0. ) ;
     cout << endl << "dedx_x   : "; for(int i=0; i<nplanes; i++)cout << dedx_x[i] << " ";
     cout << endl << "dedx_y   : "; for(int i=0; i<nplanes; i++)cout << dedx_y[i] << " ";
 
+    cout << endl << "=================";
     cout << endl;
 }
 /**
  * Reset the fit parameters 
  */
-void ExtTrack::FitReset(){
+void ExtTrack::ResetFit(){
+
+  // cout << "  ExtTrack::ResetFit() "<<endl;
     for(int i=0; i<5; i++) al[i]=-9999.;
     chi2=0.;
     nstep=0;
-    for(int i=0; i<nplanes; i++) xv[i]=0.;
-    for(int i=0; i<nplanes; i++) yv[i]=0.;
-    for(int i=0; i<nplanes; i++) zv[i]=0.;
-    for(int i=0; i<nplanes; i++) axv[i]=0.;
-    for(int i=0; i<nplanes; i++) ayv[i]=0.;
+
     for(int i=0; i<5; i++) {
 	for(int j=0; j<5; j++) coval[i][j]=0.;
     }
+
+    for(int i=0; i<nplanes; i++){
+      xv[i] = 0.;
+      yv[i] = 0.;
+      zv[i] = 0.;
+      axv[i] = 0.;
+      ayv[i] = 0.;
+    }
+
+    int ngf = TrkParams::nGF;
+    for(int i=0; i<ngf; i++){
+      xGF[i] = 0.;
+      yGF[i] = 0.;
+    }
+
+
 }
 /**
- * Method to fill minimization-routine common
- * NB (zini=23.5 is implicitelly set)
+ * Reset x/y measured coordinates (NB it does not reset z coordinates) 
  */
-void ExtTrack::FillMiniStruct(cMini2track& track){
-
-    for(int i=0; i<nplanes; i++){ 
+void ExtTrack::ResetXY(){
 
-	track.xgood[i]=XGood(i);
-	track.ygood[i]=YGood(i);
-	
-	track.xm[i]=xm[i];
-	track.ym[i]=ym[i];
-	track.zm[i]=zm[i];
-	
-	float segment = 100.; //temporaneo
-	track.xm_a[i]=xm[i];
-	track.xm_b[i]=xm[i];
-	track.ym_a[i]=ym[i];
-	track.ym_b[i]=ym[i];
-	if(       XGood(i) && !YGood(i) ){
-	    track.ym_a[i] = track.ym_a[i]+segment;
-	    track.ym_b[i] = track.ym_b[i]-segment;
-	}else if( !XGood(i) && YGood(i)){
-	    track.xm_a[i] = track.xm_a[i]+segment;
-	    track.xm_b[i] = track.xm_b[i]-segment;
-	}
+  for(int i=0; i<nplanes; i++)ResetXY(i);
+      
+    
+}
+/**
+ * \brief Tracking method. It calls F77 mini routine.
+ *
+ * @see void TrkTrack::Fit(double, int&, int, int)
+ */
+void ExtTrack::Fit(double pfixed, int& fail, int iprint){
 
-	
-       	track.resx[i]=resx[i];
-	track.resy[i]=resy[i];
-    }
+  TrkParams::Load(1); 
+  if( !TrkParams::IsLoaded(1) ){
+    cout << "void ExtTrack::Fit(double,int&,int,int) ---ERROR--- m.field not loaded "<<endl;
+    return;
+  }
+ 
+  float al_ini[] = {0.,0.,0.,0.,0.};
+  
+  fail = 0;
+  
+  FillMiniStruct();   //fill F77 common
+    
 
-    for(int i=0; i<5; i++) track.al[i]=al[i]; 
+  // if fit variables have been reset, evaluate the initial guess
+  if(al[0]==-9999.&&al[1]==-9999.&&al[2]==-9999.&&al[3]==-9999.&&al[4]==-9999.){
+    if(iprint==2)cout << "calling guessext_()"<<endl;
+    guessext_();
+  }
+  // --------------------- free momentum 
+  if(pfixed==0.) {
+    exttrack_.pfixed=0.; 
+  }
+  // --------------------- fixed momentum 
+  if(pfixed!=0.) {
+    al[4]=1./pfixed;     
+    exttrack_.pfixed=pfixed;
+  }
+  
+  //  store temporarily the initial guess
+  for(int i=0; i<5; i++) al_ini[i]=exttrack_.al[i]; 
+
+  if(iprint==2){
+    cout<<endl<<" init. guess: ";
+    for(int i=0; i<5; i++)cout << al_ini[i]<<" ";
+    cout<<endl;
+  }
+  //  ------------------------------------------
+  //  call mini routine
+  //  ------------------------------------------
+  int istep=0;
+  int ifail=0;
+  if(iprint==2)cout << "calling miniext_("<<istep<<","<<ifail<<", "<<iprint<<")"<<endl;
+  miniext_(&istep,&ifail, &iprint);
+  if(ifail!=0) {
+    if(iprint)cout << "ERROR: ifail= " << ifail << endl;
+    fail = 1;
+  }
+  if(chi2!=chi2){
+    if(iprint)cout << "ERROR: chi2= " << chi2 << endl;	
+    ResetFit();
+    fail = 1;	
+  }
+  //  ------------------------------------------
+  
+  SetFromMiniStruct(); //copy content of F77 common into the ExtTrack object
+  
+  if(fail){
+    if(iprint)cout << " >>>> fit failed "<<endl;
+    for(int i=0; i<5; i++) al[i]=al_ini[i]; 
+  }
+};
 
-    track.zini = 55.;
+/**
+ * Method to fill minimization-routine common
+ * NB (zini=23.5 is implicitelly set)
+ */
+void ExtTrack::FillMiniStruct(cMiniExtTrack& track){
 
-//    track.zini = 23.5; 
-// ZINI = 23.5 !!! it should be the same parameter in all codes
+  track.nplanes = nplanes;
+  for(int i=0; i<nplanes; i++){ 
+    
+    track.xgood[i] = (double) XGood(i);
+    track.ygood[i] = (double) YGood(i);
+    
+    track.xm[i] = (double) xm[i];
+    track.ym[i] = (double) ym[i];
+    track.zm[i] = (double) zm[i];
+    
+    
+    track.xm_a[i] = (double) xma[i];
+    track.xm_b[i] = (double) xmb[i];
+    track.ym_a[i] = (double) yma[i];
+    track.ym_b[i] = (double) ymb[i];
+    track.zm_a[i] = (double) zma[i];
+    track.zm_b[i] = (double) zmb[i];
+    
+    track.resx[i] = (double) resx[i];
+    track.resy[i] = (double) resy[i];
+  }
+  
+  for(int i=0; i<5; i++) track.al[i] = (double) al[i]; 
+  
+  track.zini = (double) zini;
+  
     
 }
 /**
  * Method to set values from  minimization-routine common
  */
-void ExtTrack::SetFromMiniStruct(cMini2track *track){
+void ExtTrack::SetFromMiniStruct(cMiniExtTrack *track){
 
     for(int i=0; i<5; i++) {
-	al[i]=track->al[i];
-	for(int j=0; j<5; j++) coval[i][j]=track->cov[i][j];
+      al[i] = (float) (track->al[i]);
+      for(int j=0; j<5; j++) coval[i][j] = (float) (track->cov[i][j]);
     }
-    chi2  = track->chi2;
-    nstep = track->nstep;
-    for(int i=0; i<6; i++){ 
-	xv[i]  = track->xv[i];
-	yv[i]  = track->yv[i];
-	zv[i]  = track->zv[i];
-	xm[i]  = track->xm[i];
-	ym[i]  = track->ym[i];
-	zm[i]  = track->zm[i];
-	axv[i] = track->axv[i];
-	ayv[i] = track->ayv[i];	
+    chi2  = (float) (track->chi2);
+    nstep = (float) (track->nstep);
+    for(int i=0; i<nplanes; i++){ 
+	xv[i]  = (float) (track->xv[i]);
+	yv[i]  = (float) (track->yv[i]);
+	zv[i]  = (float) (track->zv[i]);
+	xm[i]  = (float) (track->xm[i]);
+	ym[i]  = (float) (track->ym[i]);
+	zm[i]  = (float) (track->zm[i]);
+	xma[i]  = (float) (track->xm_a[i]);
+	yma[i]  = (float) (track->ym_a[i]);
+	zma[i]  = (float) (track->zm_a[i]);
+	xmb[i]  = (float) (track->xm_b[i]);
+	ymb[i]  = (float) (track->ym_b[i]);
+	zmb[i]  = (float) (track->zm_b[i]);
+	axv[i] = (float) (track->axv[i]);
+	ayv[i] = (float) (track->ayv[i]);	
     }
-    
+    zini = (float) (track->zini);
+}
+/**
+ * Method to retrieve ID (0,1,...) of x-cluster (if any) associated to this track.
+ * If no cluster is associated, ID=-1.
+ *
+ */
+int ExtTrack::GetClusterX_ID(int ip){
+    return ((int)fabs(xgood[ip]))%10000000-1;
+};
+/**
+ * Method to retrieve ID (0-xxx) of y-cluster (if any) associated to this track.
+ * 
+ */
+int ExtTrack::GetClusterY_ID(int ip){
+    return ((int)fabs(ygood[ip]))%10000000-1;
+};
+
+/**
+ * Method to retrieve the dE/dx measured on a tracker view.
+ * @param ip plane (0-5)
+ * @param iv view (0=x 1=y)
+ */
+Float_t ExtTrack::GetDEDX(int ip, int iv){
+    if(iv==0 && ip>=0 && ip<6)return fabs(dedx_x[ip]);
+    else if(iv==1 && ip>=0 && ip<6)return fabs(dedx_y[ip]);
+    else {
+	cout << "TrkTrack::GetDEDX(int ip, int iv) -- wrong input parameters "<<ip<<iv<<endl;
+	return 0.;
+    }
+}
+/**
+ * Method to evaluate the dE/dx measured on a tracker plane. 
+ * The two measurements on x- and y-view are averaged.
+ * @param ip plane (0-5)
+ */
+Float_t ExtTrack::GetDEDX(int ip){
+    if( (Int_t)XGood(ip)+(Int_t)YGood(ip) == 0 ) return 0; 
+    return (GetDEDX(ip,0)+GetDEDX(ip,1))/((Int_t)XGood(ip)+(Int_t)YGood(ip));
+};
+
+/**
+ * Method to evaluate the dE/dx averaged over all planes. 
+ */
+Float_t ExtTrack::GetDEDX(){
+    Float_t dedx=0; 
+    for(Int_t ip=0; ip<6; ip++)dedx+=GetDEDX(ip,0)*XGood(ip)+GetDEDX(ip,1)*YGood(ip); 
+    dedx = dedx/(GetNX()+GetNY());
+    return dedx;
+};
+
+/**
+ * Method to retrieve the ladder, if assigned
+ */
+int ExtTrack::GetLadder(int ip){
+    if(XGood(ip))return (int)fabs(xgood[ip]/100000000)-1;
+    if(YGood(ip))return (int)fabs(ygood[ip]/100000000)-1;
+    return -1;
+};
+/**
+ * Method to retrieve the sensor, if assigned 
+ */
+int ExtTrack::GetSensor(int ip){
+    if(XGood(ip))return (int)((int)fabs(xgood[ip]/10000000)%10)-1;
+    if(YGood(ip))return (int)((int)fabs(ygood[ip]/10000000)%10)-1;
+    return -1;
+};
+/**
+ * 
+ */
+Int_t ExtTrack::GetClusterX_Multiplicity(int ip){ 
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::GetClusterX_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return -1;
+  }
+  return (Int_t)(multmaxx[ip]/10000); 
+}
+/**
+ * 
+ */
+
+Int_t ExtTrack::GetClusterY_Multiplicity(int ip){ 
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::GetClusterY_Multiplicity("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return -1;
+  }
+  return (Int_t)(multmaxy[ip]/10000); 
+}
+/**
+ * 
+ */
+
+Int_t ExtTrack::GetClusterX_MaxStrip(int ip){ 
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::GetClusterX_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return -1;
+  }
+  return (Int_t)(multmaxx[ip]%10000); 
+}
+/**
+ * 
+ */
+
+Int_t ExtTrack::GetClusterY_MaxStrip(int ip){ 
+  if(ip<0 || ip>=nplanes){
+    cout << " ExtTrack::GetClusterY_MaxStrip("<<ip<<",...) ---> exceed vector dimension = "<<nplanes<<endl;
+    return -1;
+  }
+  return (Int_t)(multmaxy[ip]%10000); 
+}
+Float_t ExtTrack::GetRigidity(){ 
+	Float_t rig=0;
+	if(chi2>0)rig=1./al[4];
+	if(rig<0) rig=-rig;
+	return rig;
+};
+//
+Float_t ExtTrack::GetDeflection(){
+	Float_t def=0;
+	if(chi2>0)def=al[4];
+	return def;
+};
+
+
+//
+// all that follows: EM porting from TrkLevel2
+//
+Bool_t ExtTrack::IsInsideAcceptance(float toll){
+    int ngf = TrkParams::nGF;
+    for(int i=0; i<ngf; i++){
+	//
+//	cout << endl << TrkParams::GF_element[i];
+	if( 
+	    TrkParams::GF_element[i].CompareTo("S11") &&
+	    TrkParams::GF_element[i].CompareTo("S12") &&
+	    TrkParams::GF_element[i].CompareTo("S21") &&
+	    TrkParams::GF_element[i].CompareTo("S22") &&
+	    TrkParams::GF_element[i].CompareTo("T1")  &&
+	    TrkParams::GF_element[i].CompareTo("CUF") &&
+	    TrkParams::GF_element[i].CompareTo("T2")  &&
+	    TrkParams::GF_element[i].CompareTo("T3")  &&
+	    TrkParams::GF_element[i].CompareTo("T4")  &&
+	    TrkParams::GF_element[i].CompareTo("T5")  &&
+	    TrkParams::GF_element[i].CompareTo("CLF") &&
+	    TrkParams::GF_element[i].CompareTo("T6")  &&
+	    TrkParams::GF_element[i].CompareTo("S31") &&
+	    TrkParams::GF_element[i].CompareTo("S32") &&
+	    true)continue;
+	// apply condition only within the cavity
+//	cout << " -- "<<xGF[i]<<" "<<yGF[i];
+	if( 
+	    xGF[i] <= TrkParams::xGF_min[i] + toll ||
+	    xGF[i] >= TrkParams::xGF_max[i] - toll ||
+	    yGF[i] <= TrkParams::yGF_min[i] + toll ||
+	    yGF[i] >= TrkParams::yGF_max[i] - toll ||
+	    false){
+	    
+	    return false;
+	}
+    }
+    return true;
+}
+
+/**
+ * Returns the reduced chi-square of track x-projection
+ */
+Float_t  ExtTrack::GetChi2X(){
+    float chiq=0;
+    for(int ip=0; ip<nplanes; ip++)if(XGood(ip))chiq+= pow((xv[ip]-xm[ip])/resx[ip],2.); 
+    if(GetNX()>3)chiq=chiq/(GetNX()-3);
+    else chiq=0;
+    if(chiq==0)cout << " Float_t  ExtTrack::GetChi2X() -- WARNING -- value not defined "<<chiq<<endl;
+    return chiq;
+}
+/**
+ * Returns the reduced chi-square of track y-projection
+ */
+Float_t  ExtTrack::GetChi2Y(){
+    float chiq=0;
+    for(int ip=0; ip<nplanes; ip++)if(YGood(ip))chiq+= pow((yv[ip]-ym[ip])/resy[ip],2.); 
+    if(GetNY()>2)chiq=chiq/(GetNY()-2);
+    else chiq=0;
+    if(chiq==0)cout << " Float_t  ExtTrack::GetChi2Y() -- WARNING -- value not defined "<<chiq<<endl;
+    return chiq;
 }
 
 ClassImp(ExtTrack);