--- DarthVader/TrackerLevel2/src/TrkLevel2.cpp	2006/05/19 13:15:54	1.1
+++ DarthVader/TrackerLevel2/src/TrkLevel2.cpp	2006/06/01 09:03:09	1.2
@@ -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*); 
 }
 //--------------------------------------
@@ -110,6 +111,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(){
 //};
 //--------------------------------------
@@ -408,18 +451,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 +513,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 +544,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);