diff --git a/tracy/tools/soltracy.cc b/tracy/tools/soltracy.cc
index 069509d01a8a4ab14a74d9b6708bce2ee714f9fb..8d8bff4b442922dd921d6edb6ff4c475bf0ca626 100644
--- a/tracy/tools/soltracy.cc
+++ b/tracy/tools/soltracy.cc
@@ -18,23 +18,19 @@ extern bool freq_map;
 //****************************************************************************************
 int main(int argc, char *argv[]) {
   
+  /* for time handling */
+  uint32_t start, stop;
+  
   /*set the random value for random generator*/
   const long seed = 1121; //the default random seed number
   iniranf(seed); //initialize the seed
   setrancut(2.0); //default value of the normal cut for the normal distribution
-
-
   // turn on globval.Cavity_on and globval.radiation to get proper synchr radiation damping
   // IDs accounted too if: wiggler model and symplectic integrator (method = 1)
   globval.H_exact = false;
-  /*The following values are set in the Read_lattice( ) in soleilcommon.cc*/
-  //  globval.quad_fringe = false;              // quadrupole fringe field
-  // globval.radiation   = false;             // synchrotron radiation
-  // globval.emittance   = false;             // emittance
-  // globval.pathlength  = false;
-  //  globval.bpm         = 0;
 
  
+  /* parameters to read the user input script .prm */
   long i=0L; //initilize the for loop to read command string
   char CommandStr[max_str];
   double nux = 0, nuy = 0, ksix = 0, ksiy = 0;
@@ -50,10 +46,6 @@ int main(int argc, char *argv[]) {
   UserCommand UserCommandFlag[500];
  
 
-  // for time handling
-  uint32_t start, stop;
-
- 
   /************************************************************************
    read in files and flags
    *************************************************************************/
@@ -65,17 +57,48 @@ int main(int argc, char *argv[]) {
     fprintf(stdout, "Not enough parameters\nSyntax is program parameter file\n");
     exit_(1);
   }
-
+     
+  /* get the family index of the special elements, prepare for printglob()*/
+  if(strcmp(bpm_name,"")==0)
+    globval.bpm = 0;
+  else
+    globval.bpm = ElemIndex(bpm_name);
+  if(strcmp(skew_quad_name,"")==0)
+    globval.qt = 0;  
+  else
+    globval.qt = ElemIndex(skew_quad_name);
+  if(strcmp(hcorr_name,"")==0)
+    globval.hcorr = 0;
+  else
+    globval.hcorr = ElemIndex(hcorr_name);
+  if(strcmp(vcorr_name,"")==0)
+    globval.vcorr = 0;
+  else
+    globval.vcorr = ElemIndex(vcorr_name);
+  if(strcmp(gs_name,"")==0)
+    globval.gs = 0;
+  else
+    globval.gs = ElemIndex(gs_name);
+  if(strcmp(ge_name,"")==0)
+    globval.ge = 0;
+  else
+    globval.ge = ElemIndex(ge_name);
+  
+//  globval.g = ElemIndex("g");  /* get family index of  girder*/
+   
+  /* print the summary of the element in lattice */  
+  printglob();   
+  
   /************************************************************************
     print files, very important file for debug
    *************************************************************************/
+  
   //print flat file with all the design values of the lattice,
   prtmfile("flat_file.dat");
-  // print close orbit to a file
+  // print location, twiss parameters and close orbit at all elements position to a file
   getcod(dP, lastpos);
   prt_cod("cod.out", globval.bpm, true);
 
-  
   //get_matching_params_scl(); // get tunes and beta functions at entrance
   // compute up to 3rd order momentum compact factor
   get_alphac2(); 
@@ -83,35 +106,25 @@ int main(int argc, char *argv[]) {
   //dnu_dA(10e-3, 5e-3, 0.002);
   //get_ksi2(0.0); // this gets the chromas and writes them into chrom2.out
   
+    
   
+   
   
-   // Flag factory
+/********************************************************************* 
+                            Flag factory 
+*********************************************************************/
+   
   for(i=0; i<=CommandNo; i++){//read user defined command by sequence
    //assign user defined command
     strcpy(CommandStr,UserCommandFlag[i].CommandStr); 
     
-    //print the twiss paramters in a file defined by the name
-    if(strcmp(CommandStr,"PrintTwissFlag") == 0) {
-      cout << "\n";
-      cout << "print the twiss parameters to file: "<< twiss_file << "\n";
-      printlatt(twiss_file);  
-   }
-   //print the close orbit
-   else if(strcmp(CommandStr,"PrintCODFlag") == 0) {
-      cout << "\n";
-      cout << "print the close orbit to file: "<< cod_file << "\n";
-      getcod(dP, lastpos);
-      prt_cod(cod_file, globval.bpm, true);
-   }
-    
-    
-    //turn on flag for quadrupole fringe field
-   else if(strcmp(CommandStr,"QuadFringeOnFlag") == 0) {
+       //turn on flag for quadrupole fringe field
+   if(strcmp(CommandStr,"QuadFringeOnFlag") == 0) {
       globval.quad_fringe = true;
       cout << "\n";
       cout << "globval.quad_fringe is " << globval.quad_fringe << "\n";
    }
-    
+   
   //set RF voltage  
   else if(strcmp(CommandStr,"RFvoltageFlag") == 0) {
     printf("\nSetting RF voltage:\n");
@@ -124,11 +137,159 @@ int main(int argc, char *argv[]) {
 
   // Chamber factory
   else if(strcmp(CommandStr,"ReadChamberFlag") == 0) {
-    ReadCh(chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
+    ReadCh(UserCommandFlag[i].chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
   PrintCh(); // print chamber into chamber.out
   }
 
- 
+   
+  
+
+  // read the misalignment errors to the elements, then do COD correction 
+   //  using SVD method.
+  //  Based on the function error_and_correction() in nsls-ii_lib_templ.h
+
+ // else if(strcmp(ReadaefileFlag, "") != 0){
+  else if (strcmp(CommandStr,"ReadaefileFlag") == 0) {
+  // *****  Apply corrections and output flatfile for n_stat sets of random #'s
+    bool    cod;
+    int     k; 
+    FILE    *hOrbitFile, *vOrbitFile ;
+    int     hcorrIdx[nCOR], vcorrIdx[nCOR]; //list of corr for orbit correction
+   
+    //initialize the corrector list 
+    for ( k = 0; k < nCOR; k++){
+    hcorrIdx[k] = -1;
+    vcorrIdx[k] = -1;
+    }
+
+    //Get response matrix between bpm and correctors, and then print the SVD setting to the files
+    if (n_orbit!=0 || n_scale!=0) {
+  
+      // select correctors to be used
+      readCorrectorList(hcorr_file, vcorr_file, hcorrIdx, vcorrIdx);
+   
+      fprintf(stdout, "\n\nSVD correction setting:\n");
+      fprintf(stdout, "H-plane %d singular values:\n", nwh);
+      fprintf(stdout, "V-plane %d singular values:\n\n",nwv);
+    
+      // compute beam response matrix
+      printf("\n");
+      printf("Computing beam response matrix\n");
+      //get the response matrix between bpm and correctors
+      gcmats(globval.bpm, globval.hcorr, 1, hcorrIdx); 
+      gcmats(globval.bpm, globval.vcorr, 2, vcorrIdx);
+      /*    gcmat(globval.bpm, globval.hcorr, 1); 
+        gcmat(globval.bpm, globval.vcorr, 2);*/
+    
+      // print response matrices to files 'svdh.out' and file 'svdv.out'
+      prt_gcmat(globval.bpm, globval.hcorr, 1);
+      prt_gcmat(globval.bpm, globval.vcorr, 2);  
+  
+    }  
+
+     //print the statistics of orbit in file 'OrbScanFile.out'
+     OrbScanFile = file_write(OrbScanFileName);  
+    
+  //write files with orbits at all element locations
+  hOrbitFile = file_write(hOrbitFileName);
+  vOrbitFile = file_write(vOrbitFileName);
+    
+  fprintf(hOrbitFile, "# First line: s-location (m) \n");
+  fprintf(hOrbitFile, "# After orbit correction:  Horizontal closed orbit at all element locations (with %3d BPMs) at different loop\n", GetnKid(globval.bpm));
+  fprintf(vOrbitFile, "# First line s-location (m) \n");
+  fprintf(vOrbitFile, "# After orbit correction:  Vertical closed orbit at all element locations (with %3d BPMs) at different loop\n", GetnKid(globval.bpm));
+   
+  for (k = 0; k < globval.Cell_nLoc; k++){
+    fprintf(hOrbitFile, "% 9.3e  ", Cell[k].S);
+    fprintf(vOrbitFile, "% 9.3e  ", Cell[k].S);
+  } // end for
+    
+  fprintf(hOrbitFile, "\n");
+  fprintf(vOrbitFile, "\n"); 
+    
+    for (k = 1; k <= n_stat; k++) {// loop for n_stat sets of errors 
+       
+      if (n_orbit!=0 || n_scale!=0) {      
+     
+        cod = 
+	CorrectCOD_Ns(hOrbitFile, vOrbitFile, UserCommandFlag[i].ae_file,
+	              n_orbit,n_scale,k,nwh,nwv, hcorrIdx, vcorrIdx);
+        /*      cod = CorrectCOD_N(ae_file, n_orbit, n_scale, k);*/
+        printf("\n");
+      	
+        if (cod){
+        /*         printf("done with orbit correction, now do coupling",
+               " correction plus vert. disp\n");*/
+          if(n_orbit == 0)
+	    printf("iter # %3d n_obit = 0, no orbit correction \n",k);
+	  else  
+	    printf("iter # %3d Orbit correction succeeded\n", k);
+        }
+	else
+          chk_cod(cod, "iter # %3d error_and_correction");
+      }
+      // Ring_GetTwiss(chroma, dp);
+      Ring_GetTwiss(true, 0.0);
+  
+      // print flat lattice
+      //       sprintf(mfile_name, "flat_file.%03d.dat",k);
+      //       prtmfile(mfile_name);
+
+    }   
+    
+    prt_cod("corr_after.out", globval.bpm, true);
+   
+    // close file giving orbit at BPM location
+    fclose(hOrbitFile);
+    fclose(vOrbitFile); 
+    fclose(OrbScanFile);
+  }
+  
+   // set the field error into the lattice 
+  // The coresponding field error is replaced by the new value.
+  // This feature is generatic, works for all lattice
+   else if (strcmp(CommandStr,"ReadfefileFlag") == 0) {
+    fprintf(stdout,"\n Read field error from fe_file: \n");
+    LoadFieldErrs(UserCommandFlag[i].fe_file, true, 1.0, true, 1);
+	
+    Ring_GetTwiss(true, 0.0);
+    prtmfile("flat_file_fefile.dat");
+  }
+  // read multipole error from a file; specfic for soleil lattice
+  else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) {
+    fprintf(stdout,"\n Read Multipoles file for lattice with thick sextupoles, specific for SOLEIL lattice: \n");
+    ReadFieldErr(multipole_file);
+    //first print the full lattice with error as a flat file
+    prtmfile("flat_file_errmultipole.dat"); // writes flat file with all element errors  /* very important file for debug*/
+    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
+    printglob();
+  }
+
+     //print the twiss paramters in a file defined by the name
+   else if(strcmp(CommandStr,"PrintTwissFlag") == 0) {
+      cout << "\n";
+      cout << "print the twiss parameters to file: "<< twiss_file << "\n";
+      printlatt(twiss_file);  
+   }
+   //print the close orbit
+   else if(strcmp(CommandStr,"PrintCODFlag") == 0) {
+      cout << "\n";
+      cout << "print the close orbit to file: "<< cod_file << "\n";
+      getcod(dP, lastpos);
+      prt_cod(cod_file, globval.bpm, true);
+   }
+    //print the girder
+  // else if(strcmp(CommandStr,"PrintGirderFlag") == 0) {
+   //   cout << "\n";
+   //   cout << "print the information of girder to file: "<< girder_file << "\n";
+      
+      
+     // getcod(dP, lastpos);
+     // prt_cod(cod_file, globval.bpm, true);
+  // }
+   
+
+  
   // compute tunes by tracking (should be the same as by tps)
   else if (strcmp(CommandStr,"TuneTracFlag") == 0) {
     GetTuneTrac(1026L, 0.0, &nux, &nuy);
@@ -282,18 +443,9 @@ int main(int argc, char *argv[]) {
     printglob(); /* print parameter list */
   }
 
-  //set multipoles in all magnets
-  // read multipole error from a file
-  else if(strcmp(CommandStr,"ReadMultipoleFlag") == 0) {
-    fprintf(stdout,"\n Read Multipoles file for lattice with thick sextupoles \n");
-    ReadFieldErr(multipole_file);
-    //first print the full lattice with error as a flat file
-    prtmfile("flat_file_errmultipole.dat"); // writes flat file with all element errors  /* very important file for debug*/
-    Ring_GetTwiss(chroma = true, 0.0); /* Compute and get Twiss parameters */
-    printglob();
-  }
 
- 
+  
+  
 
   /******************************************************************************************/
   // COMPUTATION PART after setting the model
@@ -515,7 +667,7 @@ int main(int argc, char *argv[]) {
     // finally, write the momentum acceptance to file "mom_aper.out".
     if (strcmp(CommandStr,"TousTrackFlag") == 0) {
       globval.Aperture_on = true;
-      ReadCh(chamber_file);
+      ReadCh(UserCommandFlag[i].chamber_file);
       //      LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
       tau = Touschek(Qb, globval.delta_RF, false, globval.eps[X_],
           globval.eps[Y_], sigma_delta, sigma_s, n_turns, true, sum_delta,
@@ -535,6 +687,7 @@ int main(int argc, char *argv[]) {
     printf("Wrong!!!!!");
  }//end of looking for user defined flag
 
+
   
   return 0;
 }//end of main() 
diff --git a/tracy/tracy/inc/read_script.h b/tracy/tracy/inc/read_script.h
index 081b8ac82677af2dc957508fab773058b3ebdb7c..3db863c84d7bea83cd605e20121fbbaab150fc2f 100644
--- a/tracy/tracy/inc/read_script.h
+++ b/tracy/tracy/inc/read_script.h
@@ -4,8 +4,19 @@
    public:
    char  CommandStr[max_str];   // bool flag name
    
-   /* parameters */
+ /***** parameters *******/
+    //RFVoltFlag
    double RFvolt;   // RF voltage
+   
+   //chamber file
+   char chamber_file[max_str];
+   
+    // misalignment error file
+   char    ae_file[max_str];
+   
+   // field error file
+   char    fe_file[max_str]; 
+   
     // FitTuneFlag ; 
    char qf[max_str],qd[max_str]; 
    double targetnux, targetnuz;
@@ -80,8 +91,6 @@
 /* random rotation coupling error*/
  err_seed=0L;  err_rms=0.0;  
 
-
-
 /* momentum acceptance */
  TrackDim[3] = '6D';
  _MomentumAccFlag_istart=1L, _MomentumAccFlag_istop=108L,
@@ -95,24 +104,38 @@
  _Phase_nturn=512L;
  _Phase_Dim[3]='4D';
  _Phase_Damping = false;
-
   
 /* fit tunes for full quadrupole*/
  targetnux = 0.0, targetnuz = 0.0;
 
 /* fit charomaticities*/
  targetksix = 0.0, targetksiz = 0.0;
+ 
+ 
 };  
  
  };
 
-// file names
+/***** file names************/
+ //files with multipole errors
  extern char fic_hcorr[max_str],fic_vcorr[max_str], fic_skew[max_str];
  extern char multipole_file[max_str];
- 
 
- extern char chamber_file[max_str];
+
  extern char twiss_file[max_str];
  extern char cod_file[max_str];
+ extern char girder_file[max_str];
+  
+ //files with the status of hcorr/vcorr status, to choose which correctors are used for orbit correction
+ extern char    hcorr_file[max_str], vcorr_file[max_str];
+// extern char    fe_file[max_str]; //the same as multipole_file[max_str]?????
+
+//COD correction
+ extern int nwh, nwv; 
+
+extern char hcorr_name[max_str], vcorr_name[max_str];
+extern char skew_quad_name[max_str], bpm_name[max_str];
+extern char gs_name[max_str], ge_name[max_str];
+
 // function
  void read_script(const char *param_file_name, bool rd_lat,long& CommNo, UserCommand UserCommandFlag[]);