--- DarthVader/TrackerLevel2/src/TrkLevel2.cpp	2006/05/19 13:15:54	1.1.1.1
+++ DarthVader/TrackerLevel2/src/TrkLevel2.cpp	2006/06/30 09:22:04	1.4
@@ -10,6 +10,7 @@
 //......................................
 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*); 
 }
 //--------------------------------------
@@ -17,7 +18,8 @@
 //
 //--------------------------------------
 TrkTrack::TrkTrack(){
-    image = 0; 
+    seqno = -1; 
+    image = -1; 
     chi2  = 0;
     for(int it1=0;it1<5;it1++){
 	al[it1] = 0;
@@ -46,6 +48,7 @@
 //
 //--------------------------------------
 TrkTrack::TrkTrack(const TrkTrack& t){
+    seqno = t.seqno; 
     image = t.image; 
     chi2  = t.chi2;
     for(int it1=0;it1<5;it1++){
@@ -110,6 +113,48 @@
 //
 //
 //--------------------------------------
+/**
+ * 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 TrkTrack::DoTrack2(Trajectory* t){
+
+    double *dxout   = new double[t->npoint];
+    double *dyout   = new double[t->npoint];
+    double *dthxout = new double[t->npoint];
+    double *dthyout = new double[t->npoint];
+    double *dtlout  = new double[t->npoint];
+    double *dzin    = new double[t->npoint];
+    double dal[5];
+
+    int ifail = 0;
+
+    for (int i=0; i<5; i++)         dal[i]  = (double)al[i];
+    for (int i=0; i<t->npoint; i++) dzin[i] = (double)t->z[i];
+
+    dotrack2_(&(t->npoint),dzin,dxout,dyout,dthxout,dthyout,dtlout,dal,&ifail);
+    
+    for (int i=0; i<t->npoint; i++){
+	t->x[i]   = (float)*dxout++;
+	t->y[i]   = (float)*dyout++;
+	t->thx[i] = (float)*dthxout++;
+	t->thy[i] = (float)*dthyout++;
+	t->tl[i]  = (float)*dtlout++;
+    }
+
+//    delete [] dxout;
+//    delete [] dyout;
+//    delete [] dzin;
+
+    return ifail;
+};
+//--------------------------------------
+//
+//
+//--------------------------------------
 //float TrkTrack::BdL(){
 //};
 //--------------------------------------
@@ -195,6 +240,9 @@
     Track    = new TClonesArray("TrkTrack");
     SingletX = new TClonesArray("TrkSinglet");
     SingletY = new TClonesArray("TrkSinglet");
+//    Track    = 0;
+//    Singlet = 0;
+//    SingletY = 0;
 }
 //--------------------------------------
 //
@@ -223,6 +271,10 @@
  * Fills a TrkLevel2 object with values from a struct cTrkLevel2 (to get data from F77 common).
  */
 void TrkLevel2::FillCommonVar(cTrkLevel2 *l2){
+    //
+//    Track    = new TClonesArray("TrkTrack");
+//    SingletX = new TClonesArray("TrkSinglet");
+//    SingletY = new TClonesArray("TrkSinglet");
 //  temporary objects:
     TrkSinglet* t_singlet = new TrkSinglet();
     TrkTrack*   t_track   = new TrkTrack();
@@ -234,7 +286,9 @@
 //  *** 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];
@@ -302,6 +356,11 @@
  * 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;
@@ -330,6 +389,52 @@
     }
     return sorted;
 }
+TClonesArray *TrkLevel2::GetTracks_NFitSorted(){
+
+    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){
+	int nfit =0;
+	float chi2ref=1000000;
+	// first loop to search maximum num. of fit points
+	for(int i=0; i<this->ntrk(); i++){
+	    if( ((TrkTrack *)t[i])->GetNtot() >= nfit && m[i]==1){
+		nfit =    ((TrkTrack *)t[i])->GetNtot();
+//		cout << "1** "<<i<< " " << nfit<<endl;
+	    }
+	}
+	//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;
+	N--;	
+	indo++;	
+    }
+    return sorted;
+}
 //--------------------------------------
 //
 //
@@ -408,18 +513,56 @@
 //
 //--------------------------------------
 /**
+ * Trajectory default constructor. 
+ * (By default is created with z-coordinates inside the tracking volume)
+  */
+Trajectory::Trajectory(){
+    npoint = 10;
+    x = new float[npoint];
+    y = new float[npoint];
+    z = new float[npoint];
+    thx = new float[npoint];
+    thy = new float[npoint];
+    tl = new float[npoint];
+    float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1);
+    for(int i=0; i<npoint; i++){
+	x[i] = 0;
+	y[i] = 0;
+	z[i] = (ZTRKUP) - i*dz;
+	thx[i] = 0;
+	thy[i] = 0;
+	tl[i] = 0;
+    }
+}
+//--------------------------------------
+//
+//
+//--------------------------------------
+/**
  * Trajectory constructor.
+ * (By default is created with z-coordinates inside the tracking volume)
  * \param n Number of points
  */
 Trajectory::Trajectory(int n){
+    if(n<=0){
+	cout << "NB! Trajectory must have at least 1 point >>> created with 10 points" << endl;
+	n=10;
+    }
     npoint = n;
     x = new float[npoint];
     y = new float[npoint];
     z = new float[npoint];
+    thx = new float[npoint];
+    thy = new float[npoint];
+    tl = new float[npoint];
+    float dz = ((ZTRKUP)-(ZTRKDW))/(npoint-1);
     for(int i=0; i<npoint; i++){
-	x[i] = 0;
+    	x[i] = 0;
 	y[i] = 0;
-	z[i] = 0;
+	z[i] = (ZTRKUP) - i*dz;
+	thx[i] = 0;
+	thy[i] = 0;
+	tl[i] = 0;
     }
 }
 //--------------------------------------
@@ -432,15 +575,26 @@
  * \param pz Pointer to float array, defining z coordinates
  */
 Trajectory::Trajectory(int n, float* zin){
-    npoint = n;
+    npoint = 10;
+    if(n>0)npoint = n;
     x = new float[npoint];
     y = new float[npoint];
     z = new float[npoint];
-    for(int i=0; i<npoint; i++){
+    thx = new float[npoint];
+    thy = new float[npoint];
+    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++;		
+    }while(zin[i-1] > zin[i] && i < npoint);
+    npoint=i;
+    if(npoint != n)cout << "NB! Trajectory created with "<<npoint<<" points"<<endl;
 }
 //--------------------------------------
 //
@@ -452,10 +606,33 @@
 void Trajectory::Dump(){
     cout <<endl<< "Trajectory ========== "<<endl;
     for (int i=0; i<npoint; i++){
-	cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] << endl;;
+	cout << i <<" >> " << x[i] <<" "<< y[i] <<" "<< z[i] ;
+	cout <<" -- " << thx[i] <<" "<< thy[i] ;
+	cout <<" -- " << tl[i] << endl;
     };
 }
+//--------------------------------------
+//
+//
+//--------------------------------------
+/**
+ * Get trajectory length between two points
+ * @param ifirst first point (default 0)
+ * @param ilast last point (default npoint)
+ */
+float Trajectory::GetLength(int ifirst, int ilast){
+    if( ifirst<0    ) ifirst = 0;
+    if( ilast>=npoint) ilast  = npoint-1;
+    float l=0;
+    for(int i=ifirst;i<=ilast;i++){
+	l=l+tl[i]; 
+    };
+    if(z[ilast] > ZINI)l=l-tl[ilast];
+    if(z[ifirst] < ZINI)   l=l-tl[ifirst];
+    
+    return l;
 
+}
 ClassImp(TrkLevel2);
 ClassImp(TrkSinglet);
 ClassImp(TrkTrack);