diff --git a/tracy/tracy/inc/soleillib.h b/tracy/tracy/inc/soleillib.h
index 846731a60c4d0d1af206343c6f5384326e080d96..3fc4b89e36f534a623b2a2eab9bbc6637c59e954 100644
--- a/tracy/tracy/inc/soleillib.h
+++ b/tracy/tracy/inc/soleillib.h
@@ -95,19 +95,20 @@ void set_RFVoltage(const int Fnum, const double V_RF);
 /* Read multipole errors from a file for soleil*/
 void ReadFieldErr(const char *FieldErrorFile);
 
-void AddFieldErrors(const char *name, const double r0,
-		    const int n, const double Bn, const double An);
+void AddFieldErrors(const char *name,const char *keywrd, const double r0,
+		    const int n, const double Bn, const double An); 
 		    
-void AddFieldValues_type(const int N, const double r0,
+void AddFieldValues_type(const int N, const char *keywrd, const double r0,
 			 const int n, const double Bn, const double An);
 
-void AddFieldValues_fam(const int Fnum, const double r0,
+void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
 			const int n, const double Bn, const double An);
 
-void Add_bnL_sys_elem(const int Fnum, const int Knum,
-		      const int n, const double bnL, const double anL);									 			 
-void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const double r0,
-			const int n, const double Bn, const double An);	 
+void Add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
+		      const int n, const double bnL, const double anL);
+		      								 			 
+void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const char *keywrd, const double r0,
+			const int n, const double Bn, const double An);
 			 
 /* fit tunes for soleil lattice, in which each quadrupole is cut into two parts*/			
 void FitTune4(long qf1,long qf2, long qd1, long qd2, double nux, double nuy);			 
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index d9ee18d9ccd723a3ab8826b2168ad6711b7df6fe..8acb94ebbfe38dbef72b668ade7ea84ecac91c78 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -4989,10 +4989,13 @@ void set_RFVoltage(const int Fnum, const double V_RF){
 
    Purpose:
        Read multipole errors from a file
+       
         The input format of the file is:
 	
-	keyWords    raduis when the field error is meausred "r0", field error order "n", 
-	            field error component "Bn", field error component "An"; "n","Bn,""An",...
+	seed   radom_number ; this set is optional, and only works for the rms error
+	
+	keyWords  sys/rms  raduis when the field error is meausred "r0", field error order "n", 
+	                   field error component "Bn", field error component "An"; "n","Bn,""An",...
               
    Input:
         
@@ -5016,10 +5019,11 @@ void set_RFVoltage(const int Fnum, const double V_RF){
 *****************************************************************************************************/
 void ReadFieldErr(const char *FieldErrorFile) 
 {  
-  char    in[max_str], name[max_str], *prm;
+  char    in[max_str], name[max_str],keywrd[max_str], *prm;
   char    *line;
   int     n = 0;    /* field error order*/
   int     LineNum = 0;
+  int     seed_val; // random seed number for the rms error 
   double  Bn = 0.0, An = 0.0, r0 = 0.0; /* field error components and radius when the field error is measured */
   /* conversion number from A to T.m for soleil*/
   double  _convHcorr = 8.14e-4,_convVcorr = 4.642e-4, _convQt = 93.83e-4;
@@ -5032,54 +5036,72 @@ void ReadFieldErr(const char *FieldErrorFile)
   printf("\n");
   /* read lines*/
   while (line=fgets(in, max_str, inf)) {
+  
   /* kill preceding whitespace generated by "table" key
         or "space" key, but leave \n 
         so we're guaranteed to have something*/
      while(*line == ' ' || *line == '\t') {
        line++;
      }    
+    
+    /* count line number for debug*/
     LineNum++;
+    
     /* check the line is whether comment line or null line*/
     if (strstr(line, "#") == NULL && strcmp(line,"\n") != 0 &&
          strcmp(line,"\r") != 0 &&strcmp(line,"\r\n") != 0) {
-        /*read and assign the key words and measure radius*/
+	 
+        
         sscanf(line, "%s", name); 
-	sscanf(line, " %*s %lf", &r0);
-	if (prt) printf("\nsetting multipole error to: %-5s r0 = %7.1le\n",name,r0);
-	// skip first two parameters
-	strtok(line, " \t");
-	strtok(NULL, " \t");
 	
-	/* read the end of line symbol '\n','\r','\r\n' at different operation system*/
-	while ((prm = strtok(NULL, " \t")) != NULL && strcmp(prm, "\n") != 0 &&
+	if (strcmp("seed", name) == 0) { // the line to set random seed
+          sscanf(line, "%*s %d", &seed_val); 
+          printf("ReadFieldErr: setting random seed to %d\n", seed_val);
+          iniranf(seed_val); 
+      } else{//line to set (n Bn An sequence)
+	
+	/*read and assign the key words and measure radius*/
+	  sscanf(line, " %*s %s %lf",keywrd, &r0);
+	  if (prt) printf("\nsetting multipole error to: %-5s r0 = %7.1le\n",name,r0);
+	  // skip first three parameters
+	  strtok(line, " \t");
+	  strtok(NULL, " \t");
+	  strtok(NULL, " \t");
+	
+	  /* read the end of line symbol '\n','\r','\r\n' at different operation system*/
+	  while ((prm = strtok(NULL, " \t")) != NULL && strcmp(prm, "\n") != 0 &&
 	       strcmp(prm, "\r") != 0 && strcmp(prm, "\r\n") != 0) {
 	       
-	  /* read and assign n Bn An*/
-	  sscanf(prm, "%d", &n);
-	  prm = get_prm(); /*move the pointer to the next block of the line, delimiter is table key */
-	  sscanf(prm, "%lf", &Bn);
-	  prm = get_prm(); 
-	  sscanf(prm, "%lf", &An);
+	    /* read and assign n Bn An*/
+	    sscanf(prm, "%d", &n);
+	    prm = get_prm(); /*move the pointer to the next block of the line, delimiter is table key */
+	    sscanf(prm, "%lf", &Bn);
+	    prm = get_prm(); 
+	    sscanf(prm, "%lf", &An);
 	  
 	  
 	  
-	  if (prt)
-	    printf(" n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
+	    if (prt)
+	      printf(" n = %2d, Bn = %9.1e, An = %9.1e\n", n, Bn, An);
 	  
-	  /* set multipole error to horizontal correctors of soleil ring*/
-	  if(strcmp("hcorr", name) == 0)
-	    AddCorrQtErr_fam(fic_hcorr,globval.hcorr,_convHcorr,r0,n,Bn,An);
-	  /* set multipole error to vertical correctors of soleil ring*/
-         else if(strcmp("vcorr", name) == 0)
-	    AddCorrQtErr_fam(fic_vcorr,globval.vcorr,_convVcorr,r0,n,Bn,An);
-	  /* set multipole error to skew quadrupoles of soleil ring*/
-	  else if(strcmp("qt", name) == 0)
-	    AddCorrQtErr_fam(fic_skew,globval.qt,_convQt,r0,n,Bn,An);   
-	  else
-	  /* set errors for other multipole*/
-	    AddFieldErrors(name, r0, n, Bn, An) ;
+	      
+	    /* set multipole error to horizontal correctors of soleil ring*/
+	    if(strcmp("hcorr", name) == 0)
+	      AddCorrQtErr_fam(fic_hcorr,globval.hcorr,_convHcorr,keywrd,r0,n,Bn,An);
+	    /* set multipole error to vertical correctors of soleil ring*/
+           else if(strcmp("vcorr", name) == 0)
+	     AddCorrQtErr_fam(fic_vcorr,globval.vcorr,_convVcorr,keywrd,r0,n,Bn,An);
+	    /* set multipole error to skew quadrupoles of soleil ring*/
+	     else if(strcmp("qt", name) == 0)
+	       AddCorrQtErr_fam(fic_skew,globval.qt,_convQt,keywrd,r0,n,Bn,An);   
+	    else
+	    /* set errors for other multipole*/
+	      AddFieldErrors(name,keywrd, r0, n, Bn, An) ;
 	}      
-    } else
+    }//end of read the (n Bn An) sequence
+  
+  //end of the line
+  }else
       continue;
      // printf("%s", line);
   }
@@ -5088,7 +5110,7 @@ void ReadFieldErr(const char *FieldErrorFile)
 }
 
 /***********************************************************************
-void AddFieldErrors(const char *name, const double r0,
+void AddFieldErrors(const char *name, const char *keywrd,const double r0,
 		    const int n, const double Bn, const double An) 
 		    
    Purpose:
@@ -5098,6 +5120,9 @@ void AddFieldErrors(const char *name, const double r0,
    
    Input:
       name         type name or element name
+      keywrd       "rms" or "sys"
+                   "rms":  random  multipole error
+		   "sys":  systematic multipole error
       r0           radius at which error is measured, error field is relative 
                    to the design field strength when r0 !=0
       n            order of the error
@@ -5120,7 +5145,7 @@ void AddFieldErrors(const char *name, const double r0,
  Comments:
      10/2010  Written by Jianfeng Zhang
 **********************************************************************/
-void AddFieldErrors(const char *name, const double r0,
+void AddFieldErrors(const char *name,const char *keywrd, const double r0,
 		    const int n, const double Bn, const double An) 
 {
   int     Fnum = 0;
@@ -5128,15 +5153,15 @@ void AddFieldErrors(const char *name, const double r0,
   if (strcmp("all", name) == 0) {
     printf("all: not yet implemented\n");
   } else if (strcmp("dip", name) == 0) {
-    AddFieldValues_type(Dip, r0, n, Bn, An);
+    AddFieldValues_type(Dip,keywrd, r0, n, Bn, An);
   } else if (strcmp("quad", name) == 0) {
-    AddFieldValues_type(Quad, r0, n, Bn, An);
+    AddFieldValues_type(Quad, keywrd,r0, n, Bn, An);
   } else if (strcmp("sext", name) == 0) {
-    AddFieldValues_type(Sext, r0, n, Bn, An);
+    AddFieldValues_type(Sext, keywrd, r0, n, Bn, An);
   } else {/*add error to elements*/
     Fnum = ElemIndex(name);
     if(Fnum > 0)
-      AddFieldValues_fam(Fnum, r0, n, Bn, An);
+      AddFieldValues_fam(Fnum,keywrd, r0, n, Bn, An);
     else 
       printf("SetFieldErrors: undefined element %s\n", name);
   }
@@ -5144,16 +5169,20 @@ void AddFieldErrors(const char *name, const double r0,
 
 
 /***********************************************************************
-void SetFieldValues_type(const int N, const double r0,
+void SetFieldValues_type(const int N, const char *keywrd, const double r0,
 			 const int n, const double Bn, const double An)
 			 
    Purpose:
        Add the field error of the upright multipole with the design order "type" 
        with the previous value, and then the summation value replaces the previous value.
    Input:
-      N            type name     
+      N            type name 
+      keywrd       "rms" or "sys"
+                   "rms":  random  multipole error
+		   "sys":  systematic multipole error    
       r0           radius at which error is measured, error field is relative 
                    to the design field strength when r0 != 0
+		   if r0 == 0, the Bn and An are absolute value. 
       n            order of the error
       Bn           relative B component of  n-th error
       An           relative A component of  n-th error
@@ -5174,18 +5203,17 @@ void SetFieldValues_type(const int N, const double r0,
      
  Comments:
      14/10/2010  Written by Jianfeng Zhang
+     
+     Only works for soleil lattice, since the Q2/Q7, QP2a,b/QP7a,b are 
+     long quadrupoles, which have different multipole errors from other
+     short quadrupoles
 **********************************************************************/
-void AddFieldValues_type(const int N, const double r0,
+void AddFieldValues_type(const int N, const char *keywrd, const double r0,
 			 const int n, const double Bn, const double An)
 {
   double  bnL = 0.0, anL = 0.0, KLN = 0.0;
     int   k = 0;
 
-    
-  if (r0 == 0.0) {
-    // input is: (b_n*L), (a_n*L) ????
-      set_bnL_sys_type(N, n, Bn, An);
-  } else {
       // find the strength for multipole
       for(k = 1; k <= globval.Cell_nLoc; k++)
       {
@@ -5200,26 +5228,35 @@ void AddFieldValues_type(const int N, const double r0,
 	    
 	    
 	  //absolute integrated multipole error strength
-	  bnL = Bn/pow(r0, n-N)*KLN; 
-          anL = An/pow(r0, n-N)*KLN;
+	  if (r0 == 0){ 
+	    bnL = Bn;
+	    anL = An;
+	  }else{
+	    bnL = Bn/pow(r0, n-N)*KLN; 
+            anL = An/pow(r0, n-N)*KLN;
+	  }
            
 	  
 	  //NOT add the multipole errors of short quadrupole to long quadrupole qp2 & qp7 of soleil ring
+	    // for the lattice with quadrupoles which are cut into two halves
           if(N == 2 && strncmp(Cell[k].Elem.PName,"qp2",3)==0)
-	    Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, 0, 0);
+	    Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum,keywrd, n, 0, 0);
 	  else if(N == 2 && strncmp(Cell[k].Elem.PName,"qp7",3)==0)
-	    Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, 0, 0);
+	    Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd, n, 0, 0);
+	    // for the lattice with full quadrupoles 
+	  else if(N == 2 && strncmp(Cell[k].Elem.PName,"q2",2)==0)
+	    Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd, n, 0, 0);
+	  else if(N == 2 && strncmp(Cell[k].Elem.PName,"q7",2)==0)
+	    Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd, n, 0, 0);
 	  else
 	  //add errors to multipoles except qp2, qp7
-	     Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, n, bnL,anL);	
-	
-	}
+	     Add_bnL_sys_elem(Cell[k].Fnum, Cell[k].Knum, keywrd,n, bnL,anL);	
       } 
   }  
   
 }
 /***********************************************************************
-void AddFieldValues_fam(const int Fnum, const double r0,
+void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
 			const int n, const double Bn, const double An)
 			 
    Purpose:
@@ -5228,7 +5265,12 @@ void AddFieldValues_fam(const int Fnum, const double r0,
    
    Input:
       Fnum            family name 
+      keywrd       "rms" or "sys"
+                   "rms":  random  multipole error
+		   "sys":  systematic multipole error   
       r0              radius at which error is measured
+                      for the case of r0 ???????
+		      
       n               order of the error
       Bn              relative B component for the n-th error
       An              relative A component for the n-th error
@@ -5250,7 +5292,7 @@ void AddFieldValues_fam(const int Fnum, const double r0,
  Comments:
      10/2010  Written by Jianfeng Zhang
 **********************************************************************/
-void AddFieldValues_fam(const int Fnum, const double r0,
+void AddFieldValues_fam(const int Fnum, const char *keywrd, const double r0,
 			const int n, const double Bn, const double An)
 {
   int     loc = 0, ElemNum = 0, N = 0, k = 0;
@@ -5259,18 +5301,7 @@ void AddFieldValues_fam(const int Fnum, const double r0,
   loc = Elem_GetPos(Fnum, 1); /*element index of first kid*/
   N = Cell[loc].Elem.M->n_design;/*design field order*/
   
-  if (r0 == 0.0) {
-     if(Cell[loc].Elem.PL != 0){
-       bnL = Cell[loc].Elem.PL*Bn;
-       anL = Cell[loc].Elem.PL*An;
-     }
-     else{
-       bnL = Bn;
-       anL = An;
-     }     
-    // input is: (b_n*L), (a_n*L) 
-      set_bnL_sys_fam(Fnum, n, bnL, anL);
-  } else {
+
        // find the integrated design field strength for multipole
 	   if (Cell[loc].Elem.M->n_design == 1)
            KLN = Cell[loc].Elem.PL*Cell[loc].Elem.M->Pirho; /* dipole angle */
@@ -5278,20 +5309,23 @@ void AddFieldValues_fam(const int Fnum, const double r0,
 	   KLN = GetKLpar(Cell[loc].Fnum, Cell[loc].Knum, N);/* other multipole*/
            
 	   /* absolute integrated field strength*/
+	   if (r0 == 0){ //?????????
+	    bnL = Bn;
+	    anL = An;
+	  }else{
 	   bnL = Bn/pow(r0, n-N)*KLN; 
            anL = An/pow(r0, n-N)*KLN;
-	   
+	  } 
 	   //add absolute multipole field error for the family
            for(k = 1; k <= GetnKid(Fnum); k++){
              ElemNum = Elem_GetPos(Fnum, k);  /*get the element index*/
-             Add_bnL_sys_elem(Cell[ElemNum].Fnum, Cell[ElemNum].Knum, n, bnL, anL);
+             Add_bnL_sys_elem(Cell[ElemNum].Fnum, Cell[ElemNum].Knum,keywrd, n, bnL, anL);
             }
-    }
 }
 
 
 /***********************************************************************
-void add_bnL_sys_elem(const int Fnum, const int Knum,
+void add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
 		      const int n, const double bnL, const double anL)
 			 
    Purpose:
@@ -5302,6 +5336,9 @@ void add_bnL_sys_elem(const int Fnum, const int Knum,
    Input:
       Fnum           family index
       Knum           kids index   
+      keywrd         "rms" or "sys"
+                     "rms":  random  multipole error
+		     "sys":  systematic multipole error 
       n              order of the error
       bnL            absolute integrated B component for the n-th error
       anL            absolute integrated A component for the n-th error
@@ -5322,30 +5359,52 @@ void add_bnL_sys_elem(const int Fnum, const int Knum,
      
  Comments:
      10/2010  Written Jianfeng Zhang
+     
+     Rms error on the quadrupoles: only works for full quadrupole, not for half quadrupole
 **********************************************************************/
-void Add_bnL_sys_elem(const int Fnum, const int Knum,
+void Add_bnL_sys_elem(const int Fnum, const int Knum, const char *keywrd,
 		      const int n, const double bnL, const double anL)
 {
   elemtype  elem;
+  double *elemMPB; //skew components of the multipole
+ // double *elemMPBb; //right components of the multipole
   const bool  prt = false;
 
   elem = Cell[Elem_GetPos(Fnum, Knum)].Elem;
+  
+   
+   if(strcmp("sys",keywrd)==0){ 
+
+     elemMPB = elem.M->PBsys;
+
+   }
+  if(strcmp("rms",keywrd)==0){
+    
+    elemMPB = elem.M->PBrms;
+   /* save the random scale factor of rms error PBrms*/ 
+    elem.M->PBrnd[HOMmax+n] = normranf(); 
+    elem.M->PBrnd[HOMmax-n] = normranf();  
+      
+  }
 
   if (elem.PL != 0.0) {
-    elem.M->PBsys[HOMmax+n] += bnL/elem.PL;
-    elem.M->PBsys[HOMmax-n] += anL/elem.PL;
+  
+    elemMPB[HOMmax+n] += bnL/elem.PL;
+    elemMPB[HOMmax-n] += anL/elem.PL;
   } else {
+ 
     // thin kick
-    elem.M->PBsys[HOMmax+n] += bnL; elem.M->PBsys[HOMmax-n] += anL;
+    elemMPB[HOMmax+n] += bnL; 
+    elemMPB[HOMmax-n] += anL;
   }
   
   Mpole_SetPB(Fnum, Knum, n);    //set for Bn component
   Mpole_SetPB(Fnum, Knum, -n);   //set for An component
 
   if (prt)
-    printf("add the n=%d component of %s, bnL = %e,  %e, anL = %e,  %e\n",
-	   n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
-	   bnL, elem.M->PB[HOMmax+n],anL, elem.M->PB[HOMmax-n]);
+    printf("add the %s error:  n=%d component of %s, bnL = %e,  %e, anL = %e,  %e\n",
+	   keywrd,n, Cell[Elem_GetPos(Fnum, Knum)].Elem.PName,
+	   bnL, elemMPB[HOMmax+n],anL, elemMPB[HOMmax-n]);
 }
 
 /***********************************************************************
@@ -5387,7 +5446,7 @@ void SetCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const
 		 
        10/2010  Written by Jianfeng Zhang
 **********************************************************************/
-void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const double r0,
+void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const char *keywrd, const double r0,
 			const int n, const double Bn, const double An)
 {
   int     i = 0, N = 0, corr_index = 0;
@@ -5464,14 +5523,14 @@ void AddCorrQtErr_fam(char const *fic, const int Fnum, const double conv, const
     
     if (r0 == 0.0) {
     // input is: (b_n*L), (a_n*L) ???
-      Add_bnL_sys_elem(Cell[corrlistThick[i]].Fnum,Cell[corrlistThick[i]].Knum, n, Bn, An);
+      Add_bnL_sys_elem(Cell[corrlistThick[i]].Fnum,Cell[corrlistThick[i]].Knum,keywrd, n, Bn, An);
     } else {
         conv_strength = corr*conv/brho; 
 	// absolute integrated error field strength
         bnL = Bn/pow(r0, n-N)*conv_strength; 
         anL = An/pow(r0, n-N)*conv_strength;
 	
-        Add_bnL_sys_elem(Cell[corrlistThick[i]].Fnum,Cell[corrlistThick[i]].Knum, n, bnL, anL);	
+        Add_bnL_sys_elem(Cell[corrlistThick[i]].Fnum,Cell[corrlistThick[i]].Knum,keywrd, n, bnL, anL);	
 	 
       }
     }