--- DarthVader/TrackerLevel2/src/TrkLevel2.cpp	2006/08/04 08:18:06	1.8
+++ DarthVader/TrackerLevel2/src/TrkLevel2.cpp	2006/11/21 14:00:40	1.20
@@ -4,6 +4,7 @@
  */
 #include <TrkLevel2.h>
 #include <iostream>
+#include <math.h>
 using namespace std;
 //......................................
 // F77 routines
@@ -11,7 +12,10 @@
 extern "C" {    
     void dotrack_(int*, double*, double*, double*, double*, int*);
     void dotrack2_(int*, double*, double*, double*, double*,double*, double*, double*,int*);
-    int  readb_(const char*); 
+//    int  readb_(const char*);
+    int  readb_();
+    void mini2_(int*,int*,int*); 
+    void guess_();
 }
 //--------------------------------------
 //
@@ -21,27 +25,31 @@
     seqno = -1; 
     image = -1; 
     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;
+	for(int it2=0;it2<5;it2++)coval[it1][it2] = 0;
     };
     for(int ip=0;ip<6;ip++){
-	xgood[ip]  = 0;
-	ygood[ip]  = 0;
-	xm[ip]     = 0;
-	ym[ip]     = 0;
-	zm[ip]     = 0;
-	resx[ip]   = 0;
-	resy[ip]   = 0;
-	xv[ip]     = 0;
-	yv[ip]     = 0;
-	zv[ip]     = 0;
-	axv[ip]    = 0;
-	ayv[ip]    = 0;
-	dedx_x[ip] = 0;
-	dedx_y[ip] = 0;
-    };    
+		xgood[ip]  = 0;
+		ygood[ip]  = 0;
+		xm[ip]     = 0;
+		ym[ip]     = 0;
+		zm[ip]     = 0;
+		resx[ip]   = 0;
+		resy[ip]   = 0;
+		xv[ip]     = 0;
+		yv[ip]     = 0;
+		zv[ip]     = 0;
+		axv[ip]    = 0;
+		ayv[ip]    = 0;
+		dedx_x[ip] = 0;
+		dedx_y[ip] = 0;
+//		clx[ip]    = 0;
+//		cly[ip]    = 0;
+	};
+	clx = new TRefArray(6,0);
+	cly = new TRefArray(6,0);
 };
 //--------------------------------------
 //
@@ -51,27 +59,32 @@
     seqno = t.seqno; 
     image = t.image; 
     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 it2=0;it2<5;it2++)coval[it1][it2] = t.coval[it1][it2];
     };
     for(int ip=0;ip<6;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];
-    };    
+		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];
+	//	clx[ip]    = 0;//<<<<pointer
+	//	cly[ip]    = 0;//<<<<pointer
+	};
+	clx = new TRefArray(*(t.clx));
+	cly = new TRefArray(*(t.cly));
+	
 };
 //--------------------------------------
 //
@@ -187,16 +200,252 @@
 //--------------------------------------
 void TrkTrack::Dump(){
     cout << endl << "========== Track " ;
+    cout << endl << "seq.  n. : "<< seqno;
+    cout << endl << "image n. : "<< image;
     cout << endl << "al       : "; for(int i=0; i<5; i++)cout << al[i] << " ";
     cout << endl << "chi^2    : "<< chi2;
+    cout << endl << "n.step   : "<< nstep;
     cout << endl << "xgood    : "; for(int i=0; i<6; i++)cout << xgood[i] ;
     cout << endl << "ygood    : "; for(int i=0; i<6; i++)cout << ygood[i] ;
     cout << endl << "xm       : "; for(int i=0; i<6; i++)cout << xm[i] << " ";
     cout << endl << "ym       : "; for(int i=0; i<6; i++)cout << ym[i] << " ";
     cout << endl << "zm       : "; for(int i=0; i<6; i++)cout << zm[i] << " ";
+    cout << endl << "xv       : "; for(int i=0; i<6; i++)cout << xv[i] << " ";
+    cout << endl << "yv       : "; for(int i=0; i<6; i++)cout << yv[i] << " ";
+    cout << endl << "zv       : "; for(int i=0; i<6; i++)cout << zv[i] << " ";
+    cout << endl << "resx     : "; for(int i=0; i<6; i++)cout << resx[i] << " ";
+    cout << endl << "resy     : "; for(int i=0; i<6; i++)cout << resy[i] << " ";
+    cout << endl << "coval    : "; for(int i=0; i<5; i++)cout << coval[0][i]<<" ";
+    cout << endl << "           "; for(int i=0; i<5; i++)cout << coval[1][i]<<" ";
+    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 << "dedx_x   : "; for(int i=0; i<6; i++)cout << dedx_x[i] << " ";
     cout << endl << "dedx_y   : "; for(int i=0; i<6; i++)cout << dedx_y[i] << " ";
+    cout << endl;
+}
+/**
+ * Set the TrkTrack position measurements 
+ */
+void TrkTrack::SetMeasure(double *xmeas, double *ymeas, double *zmeas){
+    for(int i=0; i<6; i++) xm[i]=*xmeas++;
+    for(int i=0; i<6; i++) ym[i]=*ymeas++;
+    for(int i=0; i<6; i++) zm[i]=*zmeas++;
+}
+/**
+ * Set the TrkTrack position resolution 
+ */
+void TrkTrack::SetResolution(double *rx, double *ry){
+    for(int i=0; i<6; i++) resx[i]=*rx++;
+    for(int i=0; i<6; i++) resy[i]=*ry++;
+}
+/**
+ * Set the TrkTrack good measurement 
+ */
+void TrkTrack::SetGood(int *xg, int *yg){
+    for(int i=0; i<6; i++) xgood[i]=*xg++;
+    for(int i=0; i<6; i++) ygood[i]=*yg++;
 }
+
+/**
+ * Load the magnetic field
+ */
+void TrkTrack::LoadField(TString path){
+    
+    strcpy(path_.path,path.Data());
+    path_.pathlen = path.Length();
+    path_.error   = 0;
+    readb_();
+
+};
+
+/**
+ * Method to fill minimization-routine common
+ */
+void TrkTrack::FillMiniStruct(cMini2track& track){
+
+    for(int i=0; i<6; i++){ 
+
+	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];
+	
+//      --- temporaneo ----------------------------
+//      andrebbe inserita la dimensione del sensore
+	float segment = 100.;
+	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;
+	}
+//      --- temporaneo ----------------------------
+	
+	track.resx[i]=resx[i];
+	track.resy[i]=resy[i];
+    }
+
+    for(int i=0; i<5; i++) track.al[i]=al[i]; 
+    track.zini = 23.5; 
+// ZINI = 23.5 !!! it should be the same parameter in all codes
+    
+}
+/**
+ * Method to set values from  minimization-routine common
+ */
+void TrkTrack::SetFromMiniStruct(cMini2track *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];
+    }
+    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];	
+    }
+
+    
+}
+/**
+ * Tracking method. It calls F77 mini routine.
+ */
+void TrkTrack::Fit(double pfixed, int& fail, int iprint){
+
+    float al_ini[] = {0.,0.,0.,0.,0.};
+
+    extern cMini2track track_;
+    fail = 0;
+    FillMiniStruct(track_);
+
+    // 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.)guess_();
+
+    // --------------------- free momentum 
+    if(pfixed==0.) {
+	track_.pfixed=0.; 
+    }
+    // --------------------- fixed momentum 
+    if(pfixed!=0.) {
+	al[4]=1./pfixed;     
+	track_.pfixed=pfixed;
+    }
+
+    //  store temporarily the initial guess
+    for(int i=0; i<5; i++) al_ini[i]=track_.al[i]; 
+
+    //  ------------------------------------------
+    //  call mini routine
+    int istep=0;
+    int ifail=0;
+    mini2_(&istep,&ifail, &iprint);
+    if(ifail!=0) {
+	if(iprint)cout << "ERROR: ifail= " << ifail << endl;
+	fail = 1;
+    }
+    //  ------------------------------------------
+    
+    SetFromMiniStruct(&track_);
+//    cout << endl << "eta ===> " << track_.al[4] << endl;
+
+//     for(int i=0; i<5; i++) al[i]=track_.al[i];
+//     chi2=track_.chi2;
+//     nstep=track_.nstep;
+//     for(int i=0; i<6; i++) xv[i]=track_.xv[i];
+//     for(int i=0; i<6; i++) yv[i]=track_.yv[i];
+//     for(int i=0; i<6; i++) zv[i]=track_.zv[i];
+//     for(int i=0; i<6; i++) axv[i]=track_.axv[i];
+//     for(int i=0; i<6; i++) ayv[i]=track_.ayv[i];
+//     for(int i=0; i<5; i++) {
+// 	for(int j=0; j<5; j++) coval[i][j]=track_.cov[i][j];
+//     }
+
+    if(fail){
+ 	if(iprint)cout << " >>>> fit failed >>>> drawing initial par"<<endl;
+ 	for(int i=0; i<5; i++) al[i]=al_ini[i]; 
+    }
+
+};
+/*
+ * Reset the fit parameters 
+ */
+void TrkTrack::FitReset(){
+    for(int i=0; i<5; i++) al[i]=-9999.;
+    chi2=0.;
+    nstep=0;
+    for(int i=0; i<6; i++) xv[i]=0.;
+    for(int i=0; i<6; i++) yv[i]=0.;
+    for(int i=0; i<6; i++) zv[i]=0.;
+    for(int i=0; i<6; i++) axv[i]=0.;
+    for(int i=0; i<6; i++) ayv[i]=0.;
+    for(int i=0; i<5; i++) {
+	for(int j=0; j<5; j++) coval[i][j]=0.;
+    }
+}
+
+//--------------------------------------
+//
+//
+//--------------------------------------
+void TrkTrack::Clear(){
+	seqno = -1; 
+	image = -1; 
+	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;
+	};
+	for(int ip=0;ip<6;ip++){
+		xgood[ip]  = 0;
+		ygood[ip]  = 0;
+		xm[ip]     = 0;
+		ym[ip]     = 0;
+		zm[ip]     = 0;
+		resx[ip]   = 0;
+		resy[ip]   = 0;
+		xv[ip]     = 0;
+		yv[ip]     = 0;
+		zv[ip]     = 0;
+		axv[ip]    = 0;
+		ayv[ip]    = 0;
+		dedx_x[ip] = 0;
+		dedx_y[ip] = 0;
+//		clx[ip]    = 0;
+//		cly[ip]    = 0;
+	};
+	clx->Clear();
+	cly->Clear();
+};
+//--------------------------------------
+//
+//
+//--------------------------------------
+void TrkTrack::Delete(){
+	Clear();
+	clx->Delete();
+	cly->Delete();
+};
+	//--------------------------------------
+//
+//
+//--------------------------------------
+
 //--------------------------------------
 //
 //
@@ -206,6 +455,7 @@
     coord[0] = 0; 
     coord[1] = 0;
     sgnl     = 0;
+	cls      = 0;
 };
 //--------------------------------------
 //
@@ -216,6 +466,8 @@
     coord[0] = s.coord[0]; 
     coord[1] = s.coord[1];
     sgnl     = s.sgnl;
+//	cls      = 0;//<<<<pointer
+	cls      = TRef(s.cls);
 };
 //--------------------------------------
 //
@@ -233,10 +485,11 @@
 //
 //--------------------------------------
 TrkLevel2::TrkLevel2(){
-    good2    = -1;
+//    good2    = -1;
     for(Int_t i=0; i<12 ; i++){
-	crc[i] = -1;
-    };
+//	crc[i] = -1;
+		good[i] = -1;
+	};
     Track    = new TClonesArray("TrkTrack");
     SingletX = new TClonesArray("TrkSinglet");
     SingletY = new TClonesArray("TrkSinglet");
@@ -249,13 +502,13 @@
 //
 //--------------------------------------
 void TrkLevel2::Dump(){
+	
     TClonesArray &t  = *Track;
     TClonesArray &sx = *SingletX;
     TClonesArray &sy = *SingletY;
-
+	//
     cout << endl << endl << "=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-";
-    cout << endl << "good2    : " << good2;
-    cout << endl << "crc      : "; for(int i=0; i<12; i++) cout << crc[i];
+    cout << endl << "good     : "; for(int i=0; i<12; i++) cout << good[i]<<" ";
     cout << endl << "ntrk()   : " << this->ntrk() ;
     cout << endl << "nclsx()  : " << this->nclsx();
     cout << endl << "nclsy()  : " << this->nclsy();
@@ -271,67 +524,164 @@
  * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
  */
 void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2){
-    //
-//    Track    = new TClonesArray("TrkTrack");
-//    SingletX = new TClonesArray("TrkSinglet");
-//    SingletY = new TClonesArray("TrkSinglet");
-//  temporary objects:
+
+	//  temporary objects:
     TrkSinglet* t_singlet = new TrkSinglet();
     TrkTrack*   t_track   = new TrkTrack();
-//  general variables
-    good2 = l2->good2;
+
+	//  **** general variables ****
+//    good2 = l2->good2;
     for(Int_t i=0; i<12 ; i++){
-	crc[i] = l2->crc[i];
-    };
-//  *** TRACKS ***
+//		crc[i] = l2->crc[i];
+		good[i] = l2->good[i];
+	};
+	//  *** TRACKS ***
     TClonesArray &t = *Track;
     for(int i=0; i<l2->ntrk; i++){
-	t_track->seqno = i;
-	t_track->image = l2->image[i]-1;
-//	cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
-	t_track->chi2  = l2->chi2_nt[i];
-	for(int it1=0;it1<5;it1++){
-	    t_track->al[it1] = l2->al_nt[i][it1];
-	    for(int it2=0;it2<5;it2++)
-		t_track->coval[it1][it2] = l2->coval[i][it2][it1];
-	};
-	for(int ip=0;ip<6;ip++){
-	    t_track->xgood[ip]  = l2->xgood_nt[i][ip];
-	    t_track->ygood[ip]  = l2->ygood_nt[i][ip];
-	    t_track->xm[ip]     = l2->xm_nt[i][ip];
-	    t_track->ym[ip]     = l2->ym_nt[i][ip];
-	    t_track->zm[ip]     = l2->zm_nt[i][ip];
-	    t_track->resx[ip]   = l2->resx_nt[i][ip];
-	    t_track->resy[ip]   = l2->resy_nt[i][ip];
-	    t_track->xv[ip]     = l2->xv_nt[i][ip];
-	    t_track->yv[ip]     = l2->yv_nt[i][ip];
-	    t_track->zv[ip]     = l2->zv_nt[i][ip];
-	    t_track->axv[ip]    = l2->axv_nt[i][ip];
-	    t_track->ayv[ip]    = l2->ayv_nt[i][ip];
-	    t_track->dedx_x[ip] = l2->dedx_x[i][ip];
-	    t_track->dedx_y[ip] = l2->dedx_y[i][ip];
-	};
-	new(t[i]) TrkTrack(*t_track);
-	t_track->Clear();
+		t_track->seqno = i;// NBNBNBNB deve sempre essere = i
+		t_track->image = l2->image[i]-1;
+	//	cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
+		t_track->chi2  = l2->chi2_nt[i];
+		t_track->nstep = l2->nstep_nt[i];
+		for(int it1=0;it1<5;it1++){
+			t_track->al[it1] = l2->al_nt[i][it1];
+			for(int it2=0;it2<5;it2++)
+			t_track->coval[it1][it2] = l2->coval[i][it2][it1];
+		};
+		for(int ip=0;ip<6;ip++){
+			t_track->xgood[ip]  = l2->xgood_nt[i][ip];
+			t_track->ygood[ip]  = l2->ygood_nt[i][ip];
+			t_track->xm[ip]     = l2->xm_nt[i][ip];
+			t_track->ym[ip]     = l2->ym_nt[i][ip];
+			t_track->zm[ip]     = l2->zm_nt[i][ip];
+			t_track->resx[ip]   = l2->resx_nt[i][ip];
+			t_track->resy[ip]   = l2->resy_nt[i][ip];
+			t_track->xv[ip]     = l2->xv_nt[i][ip];
+			t_track->yv[ip]     = l2->yv_nt[i][ip];
+			t_track->zv[ip]     = l2->zv_nt[i][ip];
+			t_track->axv[ip]    = l2->axv_nt[i][ip];
+			t_track->ayv[ip]    = l2->ayv_nt[i][ip];
+			t_track->dedx_x[ip] = l2->dedx_x[i][ip];
+			t_track->dedx_y[ip] = l2->dedx_y[i][ip];
+//			t_track->clx[ip] = 0;
+//			t_track->cly[ip] = 0;
+		};
+		new(t[i]) TrkTrack(*t_track);
+		t_track->Clear();
     };
 //  *** SINGLETS ***
     TClonesArray &sx = *SingletX;
     for(int i=0; i<l2->nclsx; i++){
-	t_singlet->plane    = l2->planex[i];
-	t_singlet->coord[0] = l2->xs[i][0];
-	t_singlet->coord[1] = l2->xs[i][1];
-	t_singlet->sgnl     = l2->signlxs[i];
-	new(sx[i]) TrkSinglet(*t_singlet);
-	t_singlet->Clear();
+		t_singlet->plane    = l2->planex[i];
+		t_singlet->coord[0] = l2->xs[i][0];
+		t_singlet->coord[1] = l2->xs[i][1];
+		t_singlet->sgnl     = l2->signlxs[i];
+//		t_singlet->cls      = 0;
+		new(sx[i]) TrkSinglet(*t_singlet);
+		t_singlet->Clear();
     }
     TClonesArray &sy = *SingletY;
     for(int i=0; i<l2->nclsy; i++){
-	t_singlet->plane    = l2->planey[i];
-	t_singlet->coord[0] = l2->ys[i][0];
-	t_singlet->coord[1] = l2->ys[i][1];
-	t_singlet->sgnl     = l2->signlys[i];
-	new(sy[i]) TrkSinglet(*t_singlet);
-	t_singlet->Clear();
+		t_singlet->plane    = l2->planey[i];
+		t_singlet->coord[0] = l2->ys[i][0];
+		t_singlet->coord[1] = l2->ys[i][1];
+		t_singlet->sgnl     = l2->signlys[i];
+//		t_singlet->cls      = 0;
+		new(sy[i]) TrkSinglet(*t_singlet);
+		t_singlet->Clear();
+	};
+	
+	delete t_track;
+	delete t_singlet;
+}
+//--------------------------------------
+//
+//
+//--------------------------------------
+/**
+ * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common). 
+ * Ref to Level1 data (clusters) is also set.
+ */
+void TrkLevel2::SetFromLevel2Struct(cTrkLevel2 *l2, TrkLevel1 *l1){
+
+//  temporary objects:
+	TrkSinglet* t_singlet = new TrkSinglet();
+	TrkTrack*   t_track   = new TrkTrack();
+// general variables
+//	good2 = l2->good2;
+	for(Int_t i=0; i<12 ; i++){
+//		crc[i] = l2->crc[i];
+		good[i] = l2->good[i];
+	};
+// *** TRACKS ***
+	TClonesArray &t = *Track;
+	for(int i=0; i<l2->ntrk; i++){
+		t_track->seqno = i;// NBNBNBNB deve sempre essere = i
+		t_track->image = l2->image[i]-1;
+//		cout << "track "<<i<<t_track->seqno << t_track->image<<endl;
+		t_track->chi2  = l2->chi2_nt[i];
+		t_track->nstep = l2->nstep_nt[i];
+		for(int it1=0;it1<5;it1++){
+			t_track->al[it1] = l2->al_nt[i][it1];
+			for(int it2=0;it2<5;it2++)
+				t_track->coval[it1][it2] = l2->coval[i][it2][it1];
+		};
+		for(int ip=0;ip<6;ip++){
+			t_track->xgood[ip]  = l2->xgood_nt[i][ip];
+			t_track->ygood[ip]  = l2->ygood_nt[i][ip];
+			t_track->xm[ip]     = l2->xm_nt[i][ip];
+			t_track->ym[ip]     = l2->ym_nt[i][ip];
+			t_track->zm[ip]     = l2->zm_nt[i][ip];
+			t_track->resx[ip]   = l2->resx_nt[i][ip];
+			t_track->resy[ip]   = l2->resy_nt[i][ip];
+			t_track->xv[ip]     = l2->xv_nt[i][ip];
+			t_track->yv[ip]     = l2->yv_nt[i][ip];
+			t_track->zv[ip]     = l2->zv_nt[i][ip];
+			t_track->axv[ip]    = l2->axv_nt[i][ip];
+			t_track->ayv[ip]    = l2->ayv_nt[i][ip];
+			t_track->dedx_x[ip] = l2->dedx_x[i][ip];
+			t_track->dedx_y[ip] = l2->dedx_y[i][ip];
+//			cout << "traccia "<<i<<"  --  "<< ip << " "<< l2->cltrx[i][ip] <<" "<< l2->cltry[i][ip] <<" "<< t_track->xgood[ip] << t_track->ygood[ip]<<endl;
+			//-----------------------------------------------------
+//			t_track->clx[ip] = l1->GetCluster(l2->cltrx[i][ip]-1);
+//			t_track->cly[ip] = l1->GetCluster(l2->cltry[i][ip]-1);
+			if(t_track->xgood[ip])t_track->clx->AddAt(l1->GetCluster(l2->cltrx[i][ip]-1),ip);
+			if(t_track->ygood[ip])t_track->cly->AddAt(l1->GetCluster(l2->cltry[i][ip]-1),ip);
+//			if(t_track->ygood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltry[i][ip] "<<l2->cltry[i][ip]<< " l1->GetCluster(l2->cltry[i][ip]-1) "<<l1->GetCluster(l2->cltry[i][ip]-1)<<endl;
+//			if(t_track->xgood[ip])cout<<" i "<<i<<" ip "<<ip<<" l2->cltrx[i][ip] "<<l2->cltrx[i][ip]<< " l1->GetCluster(l2->cltrx[i][ip]-1) "<<l1->GetCluster(l2->cltrx[i][ip]-1)<<endl;
+			//-----------------------------------------------------
+		};
+		new(t[i]) TrkTrack(*t_track);
+		t_track->Clear();
+	};
+// *** SINGLETS ***
+	TClonesArray &sx = *SingletX;
+	for(int i=0; i<l2->nclsx; i++){
+		t_singlet->plane    = l2->planex[i];
+		t_singlet->coord[0] = l2->xs[i][0];
+		t_singlet->coord[1] = l2->xs[i][1];
+		t_singlet->sgnl     = l2->signlxs[i];
+		//-----------------------------------------------------
+//		cout << "singolo x "<<i<<"  --  "<<  l2->clsx[i] <<endl;
+		t_singlet->cls      = l1->GetCluster(l2->clsx[i]-1);
+//		cout<<" i "<<i<<" l2->clsx[i] "<<l2->clsx[i]<< " l1->GetCluster(l2->clsx[i]-1) "<<l1->GetCluster(l2->clsx[i]-1)<<endl;		
+		//-----------------------------------------------------
+		new(sx[i]) TrkSinglet(*t_singlet);
+		t_singlet->Clear();
+	}
+	TClonesArray &sy = *SingletY;
+	for(int i=0; i<l2->nclsy; i++){
+		t_singlet->plane    = l2->planey[i];
+		t_singlet->coord[0] = l2->ys[i][0];
+		t_singlet->coord[1] = l2->ys[i][1];
+		t_singlet->sgnl     = l2->signlys[i];
+		//-----------------------------------------------------
+//		cout << "singolo y "<<i<<"  --  "<<  l2->clsy[i] <<endl;
+		t_singlet->cls      = l1->GetCluster(l2->clsy[i]-1);
+//		cout<<" i "<<i<<" l2->clsy[i] "<<l2->clsy[i]<< " l1->GetCluster(l2->clsy[i]-1) "<<l1->GetCluster(l2->clsy[i]-1)<<endl;		
+		//-----------------------------------------------------
+		new(sy[i]) TrkSinglet(*t_singlet);
+		t_singlet->Clear();
 	};
 	
 	delete t_track;
@@ -344,9 +694,10 @@
 void TrkLevel2::GetLevel2Struct(cTrkLevel2 *l2) const {
   
 //  general variables
-    l2->good2 = good2 ;
+//    l2->good2 = good2 ;
     for(Int_t i=0; i<12 ; i++){
-	l2->crc[i] = crc[i];
+//	l2->crc[i] = crc[i];
+		l2->good[i] = good[i];
     };
 //  *** TRACKS ***
 
@@ -354,7 +705,8 @@
     for(Int_t i=0;i<l2->ntrk;i++){
       l2->image[i] = 1 + ((TrkTrack *)Track->At(i))->image;
       l2->chi2_nt[i] =  ((TrkTrack *)Track->At(i))->chi2;
-      for(int it1=0;it1<5;it1++){
+	  l2->nstep_nt[i] =  ((TrkTrack *)Track->At(i))->nstep;
+	  for(int it1=0;it1<5;it1++){
 	l2->al_nt[i][it1] = ((TrkTrack *)Track->At(i))->al[it1];
 	for(int it2=0;it2<5;it2++)
 	  l2->coval[i][it2][it1] = ((TrkTrack *)Track->At(i))->coval[it1][it2];
@@ -398,10 +750,11 @@
 //
 //--------------------------------------
 void TrkLevel2::Clear(){
-    good2    = -1;
+//    good2    = -1;
     for(Int_t i=0; i<12 ; i++){
-	crc[i] = -1;
-    };
+//	crc[i] = -1;
+		good[i] = -1;
+	};
 /*    Track->RemoveAll();
     SingletX->RemoveAll();
     SingletY->RemoveAll();*/
@@ -414,101 +767,23 @@
 //
 //
 //--------------------------------------
+void TrkLevel2::Delete(){
+	
+	Clear();
+	Track->Delete();
+	SingletX->Delete();
+	SingletY->Delete();
+}
+//--------------------------------------
+//
+//
+//--------------------------------------
 /**
  * Sort physical tracks and stores them in a TObjectArray, ordering by increasing chi**2 value (in case of track image, it selects the one with lower chi**2). The total number of physical tracks is given by GetNTracks() and the it-th physical track can be retrieved by means of the method GetTrack(int it).
  * This method is overridden by PamLevel2::GetTracks(), where calorimeter and TOF information is used.
  */
-// TClonesArray *TrkLevel2::GetTracks(){
-// 	TClonesArray *sorted = GetTracks_NFitSorted();
-// 	return sorted;
-// };
-
-/*TClonesArray *TrkLevel2::GetTracks_Chi2Sorted(){
-
-    TClonesArray *sorted = new TClonesArray("TrkTrack");
-    TClonesArray &t = *Track;
-    TClonesArray &ts = *sorted;
-    int N=this->ntrk();
-    vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;
-
-    int indo=0;
-    int indi=0;
-    while(N != 0){
-	float chi2ref=1000000;
-	for(int i=0; i<this->ntrk(); i++){
-	    if(((TrkTrack *)t[i])->chi2 < chi2ref && m[i]==1){
-		chi2ref = ((TrkTrack *)t[i])->chi2;
-		indi = i;
-	    }
-	}
-	if( ((TrkTrack *)t[indi])->image != -1 ){
-	    m[((TrkTrack *)t[indi])->image] = 0;
-	    N--;
-	}
-	new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]);
-	m[indi] = 0;
-	N--;	
-	indo++;	
-    }
-    return sorted;
-}*/
-/*TClonesArray *TrkLevel2::GetTracks_NFitSorted(){
-
-	TClonesArray *sorted = new TClonesArray("TrkTrack");	
-    TClonesArray &t  = *Track;
-	TClonesArray &ts = *sorted;
-//    TClonesArray &ts = *PhysicalTrack;
-	int N = ntrk();
-    vector<int> m(N); for(int i=0; i<N; i++)m[i]=1;
-//	int m[50]; for(int i=0; i<N; i++)m[i]=1;
-	
-    int indo=0;
-    int indi=0;
-    while(N != 0){
-		int nfit =0;
-		float chi2ref = numeric_limits<float>::max();
-		// first loop to search maximum num. of fit points
-		for(int i=0; i < ntrk(); i++){
-//			if(N==ntrk())cout << "** "<<i<< " " << ((TrkTrack *)t[i])->GetImageSeqNo()<< " " <<((TrkTrack *)t[i])->GetNtot() << " " <<((TrkTrack *)t[i])->chi2 << endl;
-			if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){
-			nfit =    ((TrkTrack *)t[i])->GetNtot();
-			}
-		}
-		//second loop to search minimum chi2 among selected
-		for(int i=0; i<this->ntrk(); i++){
-			if(    ((TrkTrack *)t[i])->chi2 < chi2ref 
-			&& ((TrkTrack *)t[i])->GetNtot()== nfit 
-			&& m[i]==1){
-			chi2ref = ((TrkTrack *)t[i])->chi2;
-			indi = i;
-	//		cout << "2** "<<i<< " " << nfit <<" "<<chi2ref<<endl;
-			}
-		}
-		if( ((TrkTrack *)t[indi])->HasImage() ){
-			m[((TrkTrack *)t[indi])->image] = 0;
-			N--;
-	
-	//	    Int_t nfiti=((TrkTrack *)t[((TrkTrack *)t[indi])->image  ])->GetNtot();
-	//	    Float_t chi2i=((TrkTrack *)t[((TrkTrack *)t[indi])->image  ])->chi2;
-			
-	//	    cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl;
-		}
-		new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]);
-		m[indi] = 0;
-//		cout << "SORTED "<< indo << " "<< indi << " "<< N << endl; 
-		N--;	
-		indo++;	
-    }
-	m.clear();
-    return sorted;
-//    return PhysicalTrack;
-}*/
-//TClonesArray *TrkLevel2::GetTracks_NFitSorted(){
 TRefArray *TrkLevel2::GetTracks_NFitSorted(){
 
-//	cout << "GetTracks_NFitSorted(it): new TRefArray()"<< endl;
-//	TClonesArray *sorted = new TClonesArray("TrkTrack");
-//	TClonesArray &ts = *sorted;
 	TRefArray *sorted = new TRefArray();
 	
 	TClonesArray &t  = *Track;
@@ -522,34 +797,30 @@
 	while(N != 0){
 		int nfit =0;
 		float chi2ref = numeric_limits<float>::max();
+		
 		// first loop to search maximum num. of fit points
 		for(int i=0; i < ntrk(); i++){
-//			if(N==ntrk())cout << "** "<<i<< " " << ((TrkTrack *)t[i])->GetImageSeqNo()<< " " <<((TrkTrack *)t[i])->GetNtot() << " " <<((TrkTrack *)t[i])->chi2 << endl;
 			if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){
 				nfit =    ((TrkTrack *)t[i])->GetNtot();
 			}
 		}
 		//second loop to search minimum chi2 among selected
 		for(int i=0; i<this->ntrk(); i++){
-			if(    ((TrkTrack *)t[i])->chi2 < chi2ref 
-							  && ((TrkTrack *)t[i])->GetNtot()== nfit 
-							  && m[i]==1){
+			Float_t chi2 = ((TrkTrack *)t[i])->chi2; 
+			if(chi2 < 0) chi2 = chi2*1000;
+			if(    chi2 < chi2ref
+				&& ((TrkTrack *)t[i])->GetNtot() == nfit 
+				&& m[i]==1){
 				chi2ref = ((TrkTrack *)t[i])->chi2;
 				indi = i;
-	//		cout << "2** "<<i<< " " << nfit <<" "<<chi2ref<<endl;
-							  }
-		}
+			};
+		};
 		if( ((TrkTrack *)t[indi])->HasImage() ){
 			m[((TrkTrack *)t[indi])->image] = 0;
 			N--;
 	
-	//	    Int_t nfiti=((TrkTrack *)t[((TrkTrack *)t[indi])->image  ])->GetNtot();
-	//	    Float_t chi2i=((TrkTrack *)t[((TrkTrack *)t[indi])->image  ])->chi2;
-			
 	//	    cout << "i** "<< ((TrkTrack *)t[indi])->image << " " << nfiti <<" "<<chi2i<<endl;
-		}
-//		new(ts[indo]) TrkTrack(*(TrkTrack*)t[indi]);
-//		cout << "GetTracks_NFitSorted(it): Add("<<indo<<")"<< endl;
+		};
 		sorted->Add( (TrkTrack*)t[indi] );	
 		
 		m[indi] = 0;
@@ -629,21 +900,7 @@
  * Retrieves the it-th "physical" track, sorted by the method GetNTracks(). 
  * @param it Track number, ranging from 0 to GetNTracks().
  */
-/*TrkTrack *TrkLevel2::GetTrack(int it){
-    
-    if(it >= this->GetNTracks()){
-	cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; 
-	cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl;
-	return 0;
-    }
-	TClonesArray* sorted = GetTracks();
-	TClonesArray &ts = *sorted;
-	
-//    TrkTrack *track = (TrkTrack*)ts[it];
-	TrkTrack *track = new TrkTrack( *(TrkTrack*)ts[it] ); //copia
-	sorted->Delete("all");////TEMPORANEO
-	return track;
-}*/
+
 TrkTrack *TrkLevel2::GetTrack(int it){
     
 	if(it >= this->GetNTracks()){
@@ -664,7 +921,7 @@
 		
 	Float_t ntot=0;
 	TClonesArray &t = *Track;
-	for(int i=0; i<ntrk(); i++) {
+	for(int i=0; i<ntrk(); i++) {	  
 		if( ((TrkTrack *)t[i])->GetImageSeqNo() == -1 ) ntot+=1.;
 		else ntot+=0.5;
 	}
@@ -679,31 +936,6 @@
  * Retrieves (if present) the image of the it-th "physical" track, sorted by the method GetNTracks(). 
  * @param it Track number, ranging from 0 to GetNTracks().
  */
-/*TrkTrack *TrkLevel2::GetTrackImage(int it){
-
-    if(it >= this->GetNTracks()){
-	cout << "** TrkLevel2 ** Track "<< it << "does not exits! " << endl; 
-	cout << "                Physical tracks GetNTracks() = "<< this->ntrk() << endl;
-	return 0;
-    }
-	
-	TClonesArray* sorted = GetTracks();
-	TClonesArray &ts = *sorted;
-	TrkTrack *track = (TrkTrack*)ts[it];
-//    TrkTrack *track = (TrkTrack*)(*(this->GetTracks()))[it];
-    
-	if(!track->HasImage()){
-	cout << "** TrkLevel2 ** Track "<< it << "does not have image! " << endl; 
-	return 0;
-    }
-    TrkTrack *image = (TrkTrack*)(*Track)[track->image];
-
-//	GetTracks()->Delete(); ////TEMPORANEO
-	sorted->Delete(); ////TEMPORANEO
-	
-    return image;
-    
-}*/
 TrkTrack *TrkLevel2::GetTrackImage(int it){
 
 	if(it >= this->GetNTracks()){
@@ -734,8 +966,13 @@
  * Loads the magnetic field. 
  * @param s Path of the magnetic-field files.
  */
-void TrkLevel2::LoadField(TString s){
-    readb_(s.Data());
+void TrkLevel2::LoadField(TString path){
+//
+    strcpy(path_.path,path.Data());
+    path_.pathlen = path.Length();
+    path_.error   = 0;
+    readb_();
+//
 };
 //--------------------------------------
 //
@@ -833,13 +1070,13 @@
     tl = new float[npoint];
     int i=0;
     do{
-	x[i] = 0;
-	y[i] = 0;
-	z[i] = zin[i];
-	thx[i] = 0;
-	thy[i] = 0;
-	tl[i] = 0;
-	i++;		
+		x[i] = 0;
+		y[i] = 0;
+		z[i] = zin[i];
+		thx[i] = 0;
+		thy[i] = 0;
+		tl[i] = 0;
+		i++;		
     }while(zin[i-1] > zin[i] && i < npoint);
     npoint=i;
     if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;
@@ -882,6 +1119,41 @@
 
 }
 
+/**
+ * Evaluates the trajectory in the apparatus associated to the track. 
+ * It integrates the equations of motion in the magnetic field. The magnetic field should be previously loaded ( by calling  TrkLevel2::LoadField() ), otherwise an error message is returned.  
+ * @param t pointer to an object of the class Trajectory, 
+ * which z coordinates should be previously initialized by calling the proper constructor ( Trajectory::Trajectory(int n, float* zin) ).
+ * @return error flag.
+ */
+int Trajectory::DoTrack2(float* al){
+
+    double *dxout   = new double[npoint];
+    double *dyout   = new double[npoint];
+    double *dthxout = new double[npoint];
+    double *dthyout = new double[npoint];
+    double *dtlout  = new double[npoint];
+    double *dzin    = new double[npoint];
+    double dal[5];
+
+    int ifail = 0;
+
+    for (int i=0; i<5; i++)      dal[i]  = (double)al[i];
+    for (int i=0; i<npoint; i++) dzin[i] = (double)z[i];
+
+    dotrack2_(&(npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
+    
+    for (int i=0; i<npoint; i++){
+	x[i]   = (float)*dxout++;
+	y[i]   = (float)*dyout++;
+	thx[i] = (float)*dthxout++;
+	thy[i] = (float)*dthyout++;
+	tl[i]  = (float)*dtlout++;
+    }
+
+    return ifail;
+};
+
 ClassImp(TrkLevel2);
 ClassImp(TrkSinglet);
 ClassImp(TrkTrack);