From 11af03baafe0c486805d593fd761d8858fc371df Mon Sep 17 00:00:00 2001
From: zhang <zhang@9a6e40ed-f3a0-4838-9b4a-bf418f78e88d>
Date: Thu, 23 Sep 2010 12:18:13 +0000
Subject: [PATCH] Mofications Sep. 2010 Incode coupling calculation from tracy
 2 into tracy 3, and also add some comments about the functions.

---
 tracy/tools/testtracy.cc       | 204 +++++++++++++++--------------
 tracy/tracy/inc/field.h        |   3 +-
 tracy/tracy/src/nsls-ii_lib.cc | 232 ++++++++++++++++++++++++++++++++-
 tracy/tracy/src/read_script.cc |  37 +++++-
 tracy/tracy/src/soleillib.cc   |  12 +-
 tracy/tracy/src/t2lat.cc       |  25 +++-
 6 files changed, 392 insertions(+), 121 deletions(-)

diff --git a/tracy/tools/testtracy.cc b/tracy/tools/testtracy.cc
index 750d7f6..4112fc5 100644
--- a/tracy/tools/testtracy.cc
+++ b/tracy/tools/testtracy.cc
@@ -5,7 +5,7 @@ int no_tps = ORDER; // arbitrary TPSA order is defined locally
 extern bool  freq_map;
 
 #include "tracy_lib.h"
-#include "time.h"
+//#include "time.h"
 #include <sys/time.h>
 
 //***************************************************************************************
@@ -106,7 +106,7 @@ prtmfile("flat_file.dat"); // writes flat file   /* very important file for debu
   else if (ChamberNoU20Flag == true)
      DefineChNoU20();  // using vacuum chamber setting but without undulator U20
   else if (ReadChamberFlag == true)
-     ReadCh(chamber_file); /* read vacuum chamber from a file "Apertures.dat" , soleil version*/
+     ReadCh(chamber_file); /* read vacuum chamber from a file "chamber_file" , soleil version*/
 //LoadApers("Apertures.dat", 1.0, 1.0);  /* read vacuum chamber definition for bnl */
   PrintCh();
  
@@ -286,93 +286,7 @@ prtmfile("flat_file.dat"); // writes flat file   /* very important file for debu
        printf("%d %d %d \n",t_sec*1000,t_usec/1000, t_msec);
      printf("the simulation time for phase space in tracy 3 is %ld milliseconds\n",t_msec);  
    
-  }
-  
-
-  
-  //*********************************************************************************
-  //---------------------------------------------------------------------------------------------------------------------------------
-  
-  // Delicated for max4 lattice. To load alignment error files and do correction
-  if (false) {
-    //prt_cod("cod.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
-    
-    
-    globval.bpm = ElemIndex("bpm_m");                   //definition for max4 lattice, bpm
-  //  globval.bpm = ElemIndex("bpm");                         
-    globval.hcorr = ElemIndex("corr_h"); globval.vcorr = ElemIndex("corr_v");    //definition for max4 lattice, correctors
-   // globval.hcorr = ElemIndex("ch"); globval.vcorr = ElemIndex("cv");
-    globval.gs = ElemIndex("GS"); globval.ge = ElemIndex("GE");   //definition for max4 lattice, girder maker
-    
-    
-    //compute response matrix (needed for OCO)
-    gcmat(globval.bpm, globval.hcorr, 1);  gcmat(globval.bpm, globval.vcorr, 2);
-     
-
-    //print response matrix (routine in lsoc.cc)
-    //prt_gcmat(globval.bpm, globval.hcorr, 1);  prt_gcmat(globval.bpm, globval.vcorr, 2);
-       
-    //gets response matrix, does svd, evaluates correction for N seed orbits
-    //get_cod_rms_scl(dx, dy, dr, n_seed)
-    //get_cod_rms_scl(100e-6, 100e-6, 1.0e-3, 100); //trim coils aren't reset when finished
-       
-
-    //for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
-    //CorrectCOD_N("/home/simon/projects/src/lattice/AlignErr.dat", 3, 1, 1);
-   
-    //for field errors check LoadFieldErr(in nsls_ii_lib.cc) and FieldErr.dat
-    //LoadFieldErr("/home/simon/projects/src/lattice/FieldErr.dat", true, 1.0, true);
-   
-    //for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
-    //LoadAlignTol("/home/simon/projects/src/lattice/AlignErr.dat", true, 1.0, true, 1);
-    //LoadAlignTol("/home/simon/projects/out/20091126/AlignErr.dat", true, 1.0, true, 1);
-    //prt_cod("cod_err.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
-   
-      
-    // delicated for max4 lattice
-    //load alignment errors and field errors, correct orbit, repeat N times, and get statistics
-    get_cod_rms_scl_new(1); //trim coils aren't reset when finished
-  
-    
-    //for aperture limitations use LoadApers (in nsls_ii_lib.cc) and Apertures.dat
-    //globval.Aperture_on = true;
-    //LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
-    
-  }
-
-
-
-//*******************************************************************************
-//-------------------------------------------------------------------------------------------------------------------------------------  
-
-  // Call nsls-ii_lib.cc
-  // tune shift with amplitude 
-  double delta = 0;
-  if (false) {
-    cout << endl;
-    cout << "computing tune shifts" << endl;
-    dnu_dA(20e-3, 10e-3, 0.0);
-    get_ksi2(delta); // this gets the chromas and writes them into chrom2.out
- // get_ksi2(5.0e-2); // this gets the chromas and writes them into chrom2.out
-  }
-  
-  if (false) {
-    //fmap(n_x, n_y, n_tr, x_max_FMA, y_max_FMA, 0.0, true, false);
-//    fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true, false); // always use -delta_FMA (+delta_FMA appears broken)
-    fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true); // always use -delta_FMA (+delta_FMA appears broken)
-  } else {
-    globval.Cavity_on = true; // this gives longitudinal motion
-    globval.radiation = false; // this adds ripple around long. ellipse (needs many turns to resolve damp.)
-    //globval.Aperture_on = true;
-    //LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
-    //get_dynap_scl(delta, 512);
-  }
-  
-  
-
-  
-  
-  
+  } 
   
   //
   // IBS & TOUSCHEK
@@ -380,9 +294,9 @@ prtmfile("flat_file.dat"); // writes flat file   /* very important file for debu
   int     k, n_turns;
   double  sigma_s, sigma_delta, tau, alpha_z, beta_z, gamma_z, eps[3];
   FILE    *outf;
-  const double  Qb   = 5e-9;
+  const double  Qb = 5e-9;   // e charge in one bunch 
   
-  if (false) {
+  if (TouschekFlag == true) {
     double  sum_delta[globval.Cell_nLoc+1][2];
     double  sum2_delta[globval.Cell_nLoc+1][2];
     
@@ -393,7 +307,6 @@ prtmfile("flat_file.dat"); // writes flat file   /* very important file for debu
       sum_delta[k][0] = 0.0; sum_delta[k][1] = 0.0;
       sum2_delta[k][0] = 0.0; sum2_delta[k][1] = 0.0;
     }
-    
     //sigma_delta = 7.70e-04;      //410:7.70e-4,  411:9.57e-4,  412:9.12e-4
     //globval.eps[X_] = 0.326e-9;  //410:3.26e-10, 411:2.63e-10, 412:2.01e-10
     globval.eps[Y_] = 0.008e-9;
@@ -401,12 +314,14 @@ prtmfile("flat_file.dat"); // writes flat file   /* very important file for debu
     //globval.eps[Z_] = sigma_delta*sigma_s;
     //globval.delta_RF given by cav voltage in lattice file
      //globval.delta_RF = 6.20e-2; //410:6.196e-2, 411:5.285e-2, 412:4.046e-2/5.786e-2
-    n_turns = 490;               //410:490(735), 411:503(755), 412:439(659)/529(794)
+   // n_turns = 490;               //410:490(735), 411:503(755), 412:439(659)/529(794)
+     n_turns = 40;               //410:490(735), 411:503(755), 412:439(659)/529(794)
     
-
+    
+    // get the twiss parameters
     alpha_z =
-      -globval.Ascr[ct_][ct_]*globval.Ascr[delta_][ct_]
-      - globval.Ascr[ct_][delta_]*globval.Ascr[delta_][delta_];
+             -globval.Ascr[ct_][ct_]*globval.Ascr[delta_][ct_]
+             -globval.Ascr[ct_][delta_]*globval.Ascr[delta_][delta_];
     beta_z = sqr(globval.Ascr[ct_][ct_]) + sqr(globval.Ascr[ct_][delta_]);
     gamma_z = (1+sqr(alpha_z))/beta_z;
     
@@ -434,8 +349,8 @@ prtmfile("flat_file.dat"); // writes flat file   /* very important file for debu
 	     sigma_delta, sigma_s);
     
     
-    // IBS
-    if (false) {       
+    // Intra Beam Scattering(IBS)
+    if (IBSFlag == true) {       
       // initialize eps_IBS with eps_SR
       for(k = 0; k < 3; k++)
 	eps[k] = globval.eps[k];
@@ -445,9 +360,16 @@ prtmfile("flat_file.dat"); // writes flat file   /* very important file for debu
     
 
     // TOUSCHEK TRACKING
-    if (false) {       
+    // Calculate Touschek lifetime
+    // with the momentum acceptance which is determined by 
+    // the RF acceptance delta_RF and the momentum aperture 
+    // at each element location which is tracked over n turns, 
+    //the vacuum chamber is read from the file "chamber_file"  
+    // finally, write the momentum acceptance to file "mom_aper.out".
+    if (TousTrackFlag == true) {       
       globval.Aperture_on = true;
-      LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
+      ReadCh(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,
@@ -463,4 +385,86 @@ prtmfile("flat_file.dat"); // writes flat file   /* very important file for debu
     }
   
   }
+  
+  
+  
+   /*********************************************************************************
+         Delicated for max4 lattice. To load alignment error files and do correction
+  ----------------------------------------------------------------------------------
+  *********************************************************************************/
+  if (false) {
+    //prt_cod("cod.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
+    
+    
+    globval.bpm = ElemIndex("bpm_m");                   //definition for max4 lattice, bpm
+  //  globval.bpm = ElemIndex("bpm");                         
+    globval.hcorr = ElemIndex("corr_h"); globval.vcorr = ElemIndex("corr_v");    //definition for max4 lattice, correctors
+   // globval.hcorr = ElemIndex("ch"); globval.vcorr = ElemIndex("cv");
+    globval.gs = ElemIndex("GS"); globval.ge = ElemIndex("GE");   //definition for max4 lattice, girder maker
+    
+    
+    //compute response matrix (needed for OCO)
+    gcmat(globval.bpm, globval.hcorr, 1);  gcmat(globval.bpm, globval.vcorr, 2);
+     
+
+    //print response matrix (routine in lsoc.cc)
+    //prt_gcmat(globval.bpm, globval.hcorr, 1);  prt_gcmat(globval.bpm, globval.vcorr, 2);
+       
+    //gets response matrix, does svd, evaluates correction for N seed orbits
+    //get_cod_rms_scl(dx, dy, dr, n_seed)
+    //get_cod_rms_scl(100e-6, 100e-6, 1.0e-3, 100); //trim coils aren't reset when finished
+       
+
+    //for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
+    //CorrectCOD_N("/home/simon/projects/src/lattice/AlignErr.dat", 3, 1, 1);
+   
+    //for field errors check LoadFieldErr(in nsls_ii_lib.cc) and FieldErr.dat
+    //LoadFieldErr("/home/simon/projects/src/lattice/FieldErr.dat", true, 1.0, true);
+   
+    //for alignments errors check LoadAlignTol (in nsls_ii_lib.cc) and AlignErr.dat
+    //LoadAlignTol("/home/simon/projects/src/lattice/AlignErr.dat", true, 1.0, true, 1);
+    //LoadAlignTol("/home/simon/projects/out/20091126/AlignErr.dat", true, 1.0, true, 1);
+    //prt_cod("cod_err.out", globval.bpm, true); //prints a specific closed orbit with corrector strengths
+   
+      
+    // delicated for max4 lattice
+    //load alignment errors and field errors, correct orbit, repeat N times, and get statistics
+    get_cod_rms_scl_new(1); //trim coils aren't reset when finished
+  
+    
+    //for aperture limitations use LoadApers (in nsls_ii_lib.cc) and Apertures.dat
+    //globval.Aperture_on = true;
+    //LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
+    
+  }
+
+
+
+/*******************************************************************************
+    Call nsls-ii_lib.cc
+-----------------------------------------------------------------------------  
+*******************************************************************************/
+  // 
+  // tune shift with amplitude 
+  double delta = 0;
+  if (false) {
+    cout << endl;
+    cout << "computing tune shifts" << endl;
+    dnu_dA(20e-3, 10e-3, 0.0);
+    get_ksi2(delta); // this gets the chromas and writes them into chrom2.out
+ // get_ksi2(5.0e-2); // this gets the chromas and writes them into chrom2.out
+  }
+  
+  if (false) {
+    //fmap(n_x, n_y, n_tr, x_max_FMA, y_max_FMA, 0.0, true, false);
+//    fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true, false); // always use -delta_FMA (+delta_FMA appears broken)
+    fmapdp(n_x, n_dp, n_tr, x_max_FMA, -delta_FMA, 1e-3, true); // always use -delta_FMA (+delta_FMA appears broken)
+  } else {
+    globval.Cavity_on = true; // this gives longitudinal motion
+    globval.radiation = false; // this adds ripple around long. ellipse (needs many turns to resolve damp.)
+    //globval.Aperture_on = true;
+    //LoadApers("/home/simon/projects/src/lattice/Apertures.dat", 1, 1);
+    //get_dynap_scl(delta, 512);
+  }
+    
 }
diff --git a/tracy/tracy/inc/field.h b/tracy/tracy/inc/field.h
index eff09fa..78e5efd 100644
--- a/tracy/tracy/inc/field.h
+++ b/tracy/tracy/inc/field.h
@@ -269,7 +269,8 @@ template<typename T>
 inline T sqr(const T &a)
 { return a*a; }
 
-
+//overloading operator functions, specific for 
+//the user-defined data type tps
 inline tps operator+(const tps &a, const tps &b)
 { return tps(a) += b; }
 
diff --git a/tracy/tracy/src/nsls-ii_lib.cc b/tracy/tracy/src/nsls-ii_lib.cc
index 2b70a18..7c75eb0 100644
--- a/tracy/tracy/src/nsls-ii_lib.cc
+++ b/tracy/tracy/src/nsls-ii_lib.cc
@@ -465,8 +465,8 @@ void GetEmittance(const int Fnum, const bool prt)
 //	 +globval.alpha_rad[Z_])/2.0);
 
   globval.U0 = 1e9*globval.dE*globval.Energy;
-  V_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Pvolt;
-  h_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Ph;
+  V_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Pvolt; //RF voltage
+  h_RF = Cell[Elem_GetPos(Fnum, 1)].Elem.C->Ph;  // RF cavity harmonic number
   phi0 = fabs(asin(globval.U0/V_RF));
   globval.delta_RF =
     sqrt(-V_RF*cos(M_PI-phi0)*(2.0-(M_PI-2.0*(M_PI-phi0))
@@ -1888,7 +1888,33 @@ double get_Wiggler_BoBrho(const int Fnum, const int Knum)
   return Cell[Elem_GetPos(Fnum, Knum)].Elem.W->BoBrhoV[0];
 }
 
+/****************************************************************************/
+/* void set_Wiggler_BoBrho(const int Fnum, const int Knum, const double BoBrhoV)
+
+   Purpose:
+       Set the integrated field strength for the specific element in the family 
+       of wiggers.       
+       
+   Input:
+       Fnum     family name
+       BoBrhoV  integrated field strength 
 
+   Output:
+      none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       Called by set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
+       
+****************************************************************************/
 void set_Wiggler_BoBrho(const int Fnum, const int Knum, const double BoBrhoV)
 {
   Cell[Elem_GetPos(Fnum, Knum)].Elem.W->BoBrhoV[0] = BoBrhoV;
@@ -1896,7 +1922,32 @@ void set_Wiggler_BoBrho(const int Fnum, const int Knum, const double BoBrhoV)
   Wiggler_SetPB(Fnum, Knum, Quad);
 }
 
+/****************************************************************************/
+/* void set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
+
+   Purpose:
+       Set the integrated field strength for the family of wiggers       
+
+   Input:
+       Fnum     family name
+       BoBrhoV  integrated field strength 
 
+   Output:
+      none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
+
+   Comments:
+       none
+       
+****************************************************************************/
 void set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
 {
   int  k;
@@ -1905,7 +1956,35 @@ void set_Wiggler_BoBrho(const int Fnum, const double BoBrhoV)
     set_Wiggler_BoBrho(Fnum, k, BoBrhoV);
 }
 
+/****************************************************************************/
+/* void set_ID_scl(const int Fnum, const int Knum, const double scl)
+
+   Purpose:
+       Set the scale factor for the spicified element of wigger, insertion devices, 
+       or fieldMap.
+       Called by set_ID_scl(const int Fnum, const double scl).        
+
+   Input:
+       Fnum  family name
+       Knum  kid index of Fname
+       scl   scaling factor
+
+   Output:
+      none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
 
+   Comments:
+       none
+       
+****************************************************************************/
 void set_ID_scl(const int Fnum, const int Knum, const double scl)
 {
   int           k;
@@ -1933,7 +2012,32 @@ void set_ID_scl(const int Fnum, const int Knum, const double scl)
   }
 }
 
+/****************************************************************************/
+/* void set_ID_scl(const int Fnum, const double scl)
+
+   Purpose:
+       Set the scale factor for the specified wigger, insertion devices, or fieldMap.        
+
+   Input:
+       Fnum  family name
+       scl   scaling factor
+
+   Output:
+      none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
 
+   Comments:
+       none
+       
+****************************************************************************/
 void set_ID_scl(const int Fnum, const double scl)
 {
   int  k;
@@ -2680,7 +2784,33 @@ void get_IDs(void)
     }
 }
 
+/****************************************************************************/
+/* void set_ID_scl(const int Fnum, const double scl)
+
+   Purpose:
+       Set the scale factor for all the wigger, insertion devices, or fieldMap.        
+
+   Input:
+       scl   scaling factor
+
+   Output:
+      none
+
+   Return:
+       none
+
+   Global variables:
+       none
+
+   Specific functions:
+       none
 
+   Comments:
+       See also:
+                 set_ID_scl(const int Fnum, const double scl)
+		 set_ID_scl(const int Fnum, const int knum, const double scl)
+       
+****************************************************************************/
 void set_IDs(const double scl)
 {
   int  k;
@@ -3058,7 +3188,8 @@ void ini_ID_corr(void)
 
   // Configuring quads (1 --> C means thin quads located in the middle of 1s)
   // Read Betas and Nus
-  get_SQ(); Nconstr = 4*Nsext + 2;
+  get_SQ(); 
+  Nconstr = 4*Nsext + 2;
 
   // Note, allocated vectors and matrices are deallocated in ID_corr
   Xsext = dvector(1, Nconstr); Xsext0 = dvector(1, Nconstr);
@@ -3451,7 +3582,43 @@ float f_int_Touschek(const float u)
   return f;
 } 
 
+/****************************************************************************/
+/* double Touschek(const double Qb, const double delta_RF,
+		const double eps_x, const double eps_y,
+		const double sigma_delta, const double sigma_s)
+
+   Purpose:
+       Calculate Touschek lifetime 
+       using the normal theoretical formula 
+       the momentum accpetance is RF acceptance delta_RF
+      
+   Parameters:   
+             Qb  electric charge per bunch
+       delta_RF  RF momentum acceptance
+          eps_x  emittance x
+	  eps_y  emittance y
+    sigma_delta  longitudinal beam size
+        sigma_s	   
+   
+   Input:
+       none
+               
+   Output:
+       
+
+   Return:
+      Touschek lifetime [second]
+
+   Global variables:
+       
+
+   Specific functions:
+      
 
+   Comments:
+       none
+
+****************************************************************************/
 double Touschek(const double Qb, const double delta_RF,
 		const double eps_x, const double eps_y,
 		const double sigma_delta, const double sigma_s)
@@ -3463,8 +3630,8 @@ double Touschek(const double Qb, const double delta_RF,
 
   printf("\n");
   printf("Qb = %4.2f nC, delta_RF = %4.2f%%"
-	 ", eps_x = %9.3e m.rad, eps_y = %9.3e m.rad\n",
-	 1e9*Qb, 1e2*delta_RF, eps_x, eps_y);
+	 ", eps_x = %9.3e nm.rad, eps_y = %9.3e nm.rad\n",
+	 1e9*Qb, 1e2*delta_RF, 1e9*eps_x, 1e9*eps_y);
   printf("sigma_delta = %8.2e, sigma_s = %4.2f mm\n",
 	 sigma_delta, 1e3*sigma_s);
 
@@ -3575,7 +3742,55 @@ void mom_aper(double &delta, double delta_RF, const long int k,
   } 
 }
 
+/****************************************************************************/
+/* double Touschek(const double Qb, const double delta_RF,const bool consistent,
+		const double eps_x, const double eps_y,
+		const double sigma_delta, double sigma_s,
+		const int n_turn, const bool aper_on,
+		double sum_delta[][2], double sum2_delta[][2])
+
+   Purpose:
+        Calculate the Touschek lifetime
+	using the normal theoretical formula 
+	the momentum acceptance is determined by the RF acceptance delta_RF 
+	and the momentum aperture at each element location 
+	which is tracked over n turns, 
+        the vacuum chamber limitation is read from a outside file  
+      
+   Parameters:   
+             Qb  electric charge per bunch
+       delta_RF  RF momentum acceptance
+     consistent
+          eps_x  emittance x
+	  eps_y  emittance y
+    sigma_delta  longitudinal beam size
+        sigma_s	
+	 n_turn  number of turns for tracking  
+	aper_on   must be TRUE
+	         to indicate vacuum chamber setting is read from a outside file	        
+ sum_delta[][2]	
+sum2_delta[][2]		 
+		  
+   
+   Input:
+       none
+               
+   Output:
+       
+
+   Return:
+       Touschek lifetime [second]
 
+   Global variables:
+       
+
+   Specific functions:
+      mom_aper
+
+   Comments:
+       none
+
+****************************************************************************/
 double Touschek(const double Qb, const double delta_RF,const bool consistent,
 		const double eps_x, const double eps_y,
 		const double sigma_delta, double sigma_s,
@@ -3633,8 +3848,11 @@ double Touschek(const double Qb, const double delta_RF,const bool consistent,
       curly_H0 = curly_H1;
     }
 
-    sum_delta[k][0] += delta_p; sum_delta[k][1] += delta_m;
-    sum2_delta[k][0] += sqr(delta_p); sum2_delta[k][1] += sqr(delta_m);
+    sum_delta[k][0] += delta_p; 
+    sum_delta[k][1] += delta_m;
+    sum2_delta[k][0] += sqr(delta_p); 
+    sum2_delta[k][1] += sqr(delta_m);
+    
     if (prt)
       printf("%4ld %6.2f %3.2lf %3.2lf\n",
 	     k, Cell[k].S, 1e2*delta_p, 1e2*delta_m);
diff --git a/tracy/tracy/src/read_script.cc b/tracy/tracy/src/read_script.cc
index f40914e..3d030e1 100644
--- a/tracy/tracy/src/read_script.cc
+++ b/tracy/tracy/src/read_script.cc
@@ -45,6 +45,7 @@ bool InducedAmplitudeFlag = false;
 bool CodeComparaisonFlag = false;
 bool EtaFlag = false;           
 bool PhaseSpaceFlag = false;
+bool TouschekFlag = false, IBSFlag = false, TousTrackFlag = false;
 
 char chamber_file[max_str];
 
@@ -110,7 +111,7 @@ void read_script(const char *param_file_name, bool rd_lat)
       /* chamber file */
       else if (strcmp("chamber_file", name) == 0){ 
         sscanf(line, "%*s %s", str);
-        sprintf(chamber_file,"%s%s", in_dir, str); /* add file directory*/
+        sprintf(chamber_file,"%s%s", in_dir, str); /* add file directory of the chamber file*/
       } 
       
       /* read in bool flags */
@@ -380,6 +381,40 @@ void read_script(const char *param_file_name, bool rd_lat)
           exit_(1);
         } 
       }
+       else if (strcmp("TouschekFlag", name) == 0){
+        sscanf(line, "%*s %s", str);
+        if(strcmp(str, "true") == 0)
+          TouschekFlag = true;
+        else if(strcmp(str, "false") == 0)
+          TouschekFlag = false;
+        else {
+          printf("set bool flag true or flase for TouschekFlag \n" );
+          exit_(1);
+        } 
+      }
+       else if (strcmp("IBSFlag", name) == 0){
+        sscanf(line, "%*s %s", str);
+        if(strcmp(str, "true") == 0)
+          IBSFlag = true;
+        else if(strcmp(str, "false") == 0)
+          IBSFlag = false;
+        else {
+          printf("set bool flag true or flase for IBSFlag \n" );
+          exit_(1);
+        } 
+      }
+       else if (strcmp("TousTrackFlag", name) == 0){
+        sscanf(line, "%*s %s", str);
+        if(strcmp(str, "true") == 0)
+          TousTrackFlag = true;
+        else if(strcmp(str, "false") == 0)
+          TousTrackFlag = false;
+        else {
+          printf("set bool flag true or flase for TousTrackFlag \n" );
+          exit_(1);
+        } 
+      }
+      
       else {
         printf("bad line in file %s, line %ld \n", full_param_file_name, LineNum);
         exit_(1);
diff --git a/tracy/tracy/src/soleillib.cc b/tracy/tracy/src/soleillib.cc
index cbb30b4..677ed46 100644
--- a/tracy/tracy/src/soleillib.cc
+++ b/tracy/tracy/src/soleillib.cc
@@ -3281,12 +3281,12 @@ void SetSkewQuad(void)
          Based on the version in tracy 2.
 	
    Input:
-       deb first element for momentum acceptance,  "debut" is beginning in French
-       fin last element for momentum acceptance,   "fin"   is end in French
+       deb   first element for momentum acceptance,"debut" is beginning in French
+       fin   last element for momentum acceptance,"fin"   is end in French
 
-       ep_min minimum energy deviation for positive momentum acceptance
-       ep_max maximum energy deviation for positive momentum acceptance
-       nstepp number of energy steps for positive momentum acceptance
+       ep_min   minimum energy deviation for positive momentum acceptance
+       ep_max   maximum energy deviation for positive momentum acceptance
+       nstepp   number of energy steps for positive momentum acceptance
 
        em_min minimum energy deviation for negative momentum acceptance
        em_max maximum energy deviation for negative momentum acceptance
@@ -4487,7 +4487,7 @@ void Phase3(long pos, double x,double px,double y, double py,double energy,
        Should be generalized in 6D
        17/07/03 use of M_PI instead of pi
        22/03/04 save status cavity/radiation and restore it at the end
-
+       28/07/10 modified for tracy 3    
 ****************************************************************************/            
 void Coupling_Edwards_Teng(void)
 {
diff --git a/tracy/tracy/src/t2lat.cc b/tracy/tracy/src/t2lat.cc
index 542c910..c4ddf17 100644
--- a/tracy/tracy/src/t2lat.cc
+++ b/tracy/tracy/src/t2lat.cc
@@ -1974,7 +1974,25 @@ static void SetDBN(struct LOC_Lat_DealElement *LINK)
   }
 }
 
+/**************************************************************************
+ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
+                          long *errpos_, long *lc_, long *nkw_, long *inum_,
+			  long emax__, long emin__,
+                          long kmax__, long nmax__, char *chin_, char *id_,
+			  char *ElementName,
+                          char *BlockName_, double *rnum_, bool *skipflag_,
+			  bool *rsvwd_,
+                          char *line_, Lat_symbol *sym_, alfa_ *key_,
+			  Lat_symbol *ksy_,
+                          Lat_symbol *sps_, struct LOC_Lattice_Read *LINK)
 
+   Purpose:
+             Read the lattice from file fi_,
+	      first, read the element parameters,
+	      then, read the lattice
+	      element type: drift, quadrupole,sextupole,bending,cavity,.....
+
+    *************************************************************************/
 static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
                             long *errpos_, long *lc_, long *nkw_, long *inum_,
 			    long emax__, long emin__,
@@ -2039,7 +2057,6 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       L1: Drift, L=0.30;
 
     *************************************************************************/
-
   case drfsym:
     getest__(P_expset(SET, 1 << ((long)comma)), "<comma> expected", &V);
     getest__(P_addset(P_expset(SET1, 0), (long)lsym), "<L> expected", &V);
@@ -2090,7 +2107,6 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       B: bending, L=0.70, T=10.0, T1:=5.0, T2:=5.0, K=-1.0, N=8, Method=2;
 
     *************************************************************************/
-
   case bndsym:  /*4*/
     getest__(P_expset(SET, 1 << ((long)comma)), "<, > expected", &V);
     GetSym__(&V);
@@ -2227,7 +2243,6 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       Q1: Quadrupole, L=0.5, K=2.213455, N=4, Method=4;
 
     **************************************************************/
-
   case qdsym:  /*4*/
     getest__(P_expset(SET, 1 << ((long)comma)), "<, > expected", &V);
     GetSym__(&V);
@@ -2329,7 +2344,6 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       SF: Sextupole, K=-10.236345;
 
     **************************************************************************/
-
   case sexsym:  /*4*/
     QL = 0.0;   /* L */
     QK = 0.0;   /* K */
@@ -2442,7 +2456,6 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       CAV: Cavity, Frequency = 499.95e6, Voltage=1.22e6, harnum=328;
 
     **************************************************************************/
-
   case cavsym:
     ClearHOMandDBN(&V);
     getest__(P_expset(SET, 1 << ((long)comma)), "<, > expected", &V);
@@ -2502,7 +2515,7 @@ static bool Lat_DealElement(FILE **fi_, FILE **fo_, long *cc_, long *ll_,
       WITH3->Pvolt = Vrf;   /* Voltage [V] */
       WITH3->Pfreq = Frf;   /* Frequency in Hz */
       WITH3->phi = QPhi*M_PI/180.0;
-      WITH3->Ph = harnum;
+      WITH3->Ph = harnum;  /* RF harmonic number */
       SetDBN(&V);
     } else {
       printf("Elem_nFamMax exceeded: %ld(%ld)\n",
-- 
GitLab