! S. Fartoukh 30 sep 2009
! Ref sample job for tracking with multipole field imperfection in nominal 
! LHC magnets using V6.503 at injection

!---|----1----|----2----|----3----|----4----|----5----|----6----|----7----|---8;
Option,  warn,info;
System,"rm -rf temp";
System,"mkdir temp";
System,"ln -fns /afs/cern.ch/eng/lhc/optics/V6.503 db5";
System,"ln -fns /afs/cern.ch/eng/lhc/optics/SLHCV2.0 slhc";
system,"ln -fns /afs/cern.ch/eng/lhc/optics/errors/0705 wise";
 
option,-echo,-info;

!---------------------------------------------------------------;
!           CALL IN SELECTION, MACRO's, SEQUENCE, BEAM, NRJ and OPTICS FILES;
!---------------------------------------------------------------;
!Call some macro facilities
call,file="slhc/toolkit/macro.madx";

call,file="db5/V6.5.thin.seq";

!Cycling w.r.t. to IP3 (mandatory to find closed orbit in collision in the presence of errors)
seqedit,sequence=lhcb1;flatten;cycle,start=IP3;endedit;
seqedit,sequence=lhcb2;flatten;cycle,start=IP3;endedit;

!********************************************

 mylhcbeam:=1 ;// LHC beam 1 (clockwise)
!mylhcbeam:=2 ;// LHC beam 2 (clockwise)

 NRJ:= 450.0  ; // injection
!NRJ:= 7000.0 ; // collision

!********************************************

tsplit= 5.0;   !new tune split
qx00  = 64.0;  !integer h tune
qy00  = 59.0;  !integer v tune

if (NRJ<4999.9999)
{ 
!Inj. optics in thin lens for IR1&5 (11 m in IR1/5)
call,file="db5/V6.5.thin.inj.str";
qx0 = 64.28;   qy0 = 59.31;  qprime = 2.0;
Beam,particle=proton,sequence=lhcb1,energy=450.0,NPART=%NPART,sige=4.5e-4;
Beam,particle=proton,sequence=lhcb2,energy=450.0,bv = -1,NPART=%NPART,sige=4.5e-4;
};

if (NRJ>5000.0000)
{
!Coll. optics in thin lens for IR1&5 (55 cm in IR1/5)
call,file="db5/V6.5.thin.coll.str";
qx0 = 64.31;   qy0 = 59.32;  qprime = 2.0;
Beam,particle=proton,sequence=lhcb1,energy=7000.0,NPART=%NPART,sige=1.1e-4;
Beam,particle=proton,sequence=lhcb2,energy=7000.0,bv = -1,NPART=%NPART,sige=1.1e-4;         !<<<<
};

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
!                        TEST the optics
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
/*
on_x1:=0;on_x5:=0;on_sep1:=0;on_sep5:=0;
on_x2:=0;on_x8:=0;on_sep2:=0;on_sep8:=0;on_lhcb:=0;on_alice:=0;
system,"rm -f twiss.b1";system,"rm -f twiss.b2";

Use, sequence=lhcb1;
select,flag=twiss,clear;
select, flag=twiss, pattern="IP1",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
select, flag=twiss, pattern="IP2",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
select, flag=twiss, pattern="IP5",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
select, flag=twiss, pattern="IP8",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
twiss,sequence=lhcb1,file=twiss.b1;system,"cat twiss.b1";

Use, sequence=lhcb2;
select,flag=twiss,clear;
select, flag=twiss, pattern="IP1",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
select, flag=twiss, pattern="IP2",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
select, flag=twiss, pattern="IP5",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
select, flag=twiss, pattern="IP8",   column=name,s,betx,bety,alfx,alfy,dx,dpx,mux,muy;
twiss, sequence=lhcb2,file=twiss.b2;system,"cat twiss.b2";

system,"rm -f twiss.b1";system,"rm -f twiss.b2";
on_x1:=1;on_x5:=1;on_sep1:=1;on_sep5:=1;
on_x2:=1;on_x8:=1;on_sep2:=1;on_sep8:=1;on_lhcb:=-1;on_alice:=1;

Use, sequence=lhcb1;
select,flag=twiss,clear;
select, flag=twiss, pattern="IP1",   column=name,s,x,y,px,py;
select, flag=twiss, pattern="IP2",   column=name,s,x,y,px,py;
select, flag=twiss, pattern="IP5",   column=name,s,x,y,px,py;
select, flag=twiss, pattern="IP8",   column=name,s,x,y,px,py;
twiss, sequence=lhcb1, file=twiss.b1;system,"cat twiss.b1";

Use, sequence=lhcb2;
select,flag=twiss,clear;
select, flag=twiss, pattern="IP1",   column=name,s,x,y,px,py;
select, flag=twiss, pattern="IP2",   column=name,s,x,y,px,py;
select, flag=twiss, pattern="IP5",   column=name,s,x,y,px,py;
select, flag=twiss, pattern="IP8",   column=name,s,x,y,px,py;
twiss, sequence=lhcb2, file=twiss.b2;system,"cat twiss.b2";
*/

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
!                        Sequence & Beam selection
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;

if(mylhcbeam==1) {use,  sequence=lhcb1;};
if(mylhcbeam==2) {use,  sequence=lhcb2;};

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;
!               Prepare nominal twiss tables 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++;

/*
 on_x1:=1; on_sep1:=1;
 on_x2:=1; on_sep2:=1; on_alice:=1;
 on_x5:=1; on_sep5:=1;
 on_x8:=1; on_sep8:=1; on_lhcb :=-1;
 twiss;
 x.ip1=table(twiss,IP1,x);px.ip1=table(twiss,IP1,px);
 y.ip1=table(twiss,IP1,y);px.ip1=table(twiss,IP1,py); //irrelevant when recycling the sequences
*/ 
 
 on_x1:=0; on_sep1:=0;
 on_x2:=0; on_sep2:=0; on_alice:=0;
 on_x5:=0; on_sep5:=0;
 on_x8:=0; on_sep8:=0; on_lhcb :=0;
 twiss, table=nominal;   // used by orbit correction
 beta.ip1=table(twiss,IP1,betx);

! print nominal optics parameter at the MB, MQS and MSS for
! b2, b3, b4, b5, a2 and a3 correction
 select, flag=twiss, clear;
 select, flag=twiss,pattern=MB\.   ,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,pattern=MQT\.14,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,pattern=MQT\.15,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,pattern=MQT\.16,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,pattern=MQT\.17,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,pattern=MQT\.18,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,pattern=MQT\.19,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,pattern=MQT\.20,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,pattern=MQT\.21,class=multipole,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,class=MQS                      ,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 select, flag=twiss,class=MSS                      ,column=name,k0L,k1L,betx,bety,dx,mux,muy;
 twiss,  file='temp/optics0_MB.mad';

! print nominal optics parameter at the D1, MQX and triplet corrector
! for triplet correction
 select, flag=twiss, clear;
 select, flag=twiss, pattern=MQX  , class=multipole, column=name,betx,bety,x,y;
 select, flag=twiss, pattern=MBX  , class=multipole, column=name,betx,bety,x,y;
 select, flag=twiss, pattern=MBRC , class=multipole, column=name,betx,bety,x,y;
 select, flag=twiss, pattern=MQSX                  , column=name,betx,bety,x,y;
 select, flag=twiss, pattern=MCSX                  , column=name,betx,bety,x,y;
 select, flag=twiss, pattern=MCTX                  , column=name,betx,bety,x,y;
 select, flag=twiss, pattern=MCOSX                 , column=name,betx,bety,x,y;
 select, flag=twiss, pattern=MCOX                  , column=name,betx,bety,x,y;
 select, flag=twiss, pattern=MCSSX                 , column=name,betx,bety,x,y;
 twiss,  file='temp/optics0_inser.mad';
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                   Align separation magnets
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

! Still missing


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!         Call error subroutines (nominal machine and new IT/D1)
!                         and error tables
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

! Error routine and measured error table for nominal LHC

call,file="db5/measured_errors/Msubroutines.madx";
call,file="db5/toolkit/Orbit_Routines.madx";
call,file="slhc/errors/SelectLHCMonCor_phaseI.madx";
readtable, file="db5/measured_errors/rotations_Q2_integral.tab";
call,file="slhc/errors/macro_error.madx";  ! some macros for error generation 

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                Switch ON/OFF some multipole                         
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

exec ON_ALL_MULT;
ON_LSF:=1;             ! to be set to 0 if run is performed on BOINC. 
 
 ON_A1s =  0 ; ON_A1r =  0 ; ON_B1s =  0 ; ON_B1r =  0;
 ON_A2s =  1 ; ON_A2r =  1 ; ON_B2s =  1 ; ON_B2r =  1;
 ON_A3s =  1 ; ON_A3r =  1 ; ON_B3s =  1 ; ON_B3r =  1;
 ON_A4s =  1 ; ON_A4r =  1 ; ON_B4s =  1 ; ON_B4r =  1;
 ON_A5s =  1 ; ON_A5r =  1 ; ON_B5s =  1 ; ON_B5r =  1;
 ON_A6s =  1 ; ON_A6r =  1 ; ON_B6s =  1 ; ON_B6r =  1;
 ON_A7s =  1 ; ON_A7r =  1 ; ON_B7s =  1 ; ON_B7r =  1;
 ON_A8s =  1 ; ON_A8r =  1 ; ON_B8s =  1 ; ON_B8r =  1;
 ON_A9s =  1 ; ON_A9r =  1 ; ON_B9s =  1 ; ON_B9r =  1;
 ON_A10s =  1; ON_A10r =  1; ON_B10s =  1; ON_B10r =  1;
 ON_A11s =  1; ON_A11r =  1; ON_B11s =  1; ON_B11r =  1;
 ON_A12s =  ON_LSF; ON_A12r = ON_LSF; ON_B12s = ON_LSF; ON_B12r =  ON_LSF;
 ON_A13s =  ON_LSF; ON_A13r = ON_LSF; ON_B13s = ON_LSF; ON_B13r =  ON_LSF;
 ON_A14s =  ON_LSF; ON_A14r = ON_LSF; ON_B14s = ON_LSF; ON_B14r =  ON_LSF;
 ON_A15s =  ON_LSF; ON_A15r = ON_LSF; ON_B15s = ON_LSF; ON_B15r =  ON_LSF; 
  

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Define the seed number (for statistical error assignment in the new IT/D1)
! and call the corresponding measured error table for nominal LHC magnets
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

!/*
!For the mask file
eoption,seed=%SEEDRAN;
if (NRJ<4999.999) {readtable, file="wise/injection_errors-emfqcs-%SEEDRAN.tfs" ;};
if (NRJ>5000.000) {readtable, file="wise/collision_errors-emfqcs-%SEEDRAN.tfs" ;};
!*/

/*
!For the madx file
eoption,seed=1;
if (NRJ<4999.999) {readtable, file="wise/injection_errors-emfqcs-1.tfs" ;};
if (NRJ>5000.000) {readtable, file="wise/collision_errors-emfqcs-1.tfs" ;};
*/

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!             apply field errors to MB magnets
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 call,   file="db5/measured_errors/Efcomp_MB.madx"  ;
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!     correct orbit distortion resulting from MB magnets
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if((ON_A1S)^2+(ON_A1r)^2+(ON_B1S)^2+(ON_B1r)^2 >0)
{ 
  exec,   initial_micado(4);
  exec,   initial_micado(4);
};

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!           apply field errors to all other magnets
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! Nominal LHC magnets 
 
 call,   file="db5/measured_errors/Efcomp_MBRB.madx";
 call,   file="db5/measured_errors/Efcomp_MBRC.madx";
 call,   file="db5/measured_errors/Efcomp_MBRS.madx";
 call,   file="db5/measured_errors/Efcomp_MBX.madx" ;
 call,   file="db5/measured_errors/Efcomp_MBW.madx" ;
 call,   file="db5/measured_errors/Efcomp_MBXW.madx";
on_B2S:=0;call,file="db5/measured_errors/Efcomp_MQ.madx"; on_B2S:=1;
 call,   file="db5/measured_errors/Efcomp_MQM.madx" ;
 call,   file="db5/measured_errors/Efcomp_MQMC.madx";
 call,   file="db5/measured_errors/Efcomp_MQML.madx";
 call,   file="db5/measured_errors/Efcomp_MQTL.madx";
 call,   file="db5/measured_errors/Efcomp_MQW.madx" ;
on_B2S:=0;call,   file="db5/measured_errors/Efcomp_MQX.madx"; on_B2S:=1;
 call,   file="db5/measured_errors/Efcomp_MQY.madx" ;
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!               correction of field errors in MB (compatible with V6.503 & SLHC)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 select, flag=error, clear;
 select, flag=error, pattern=MB\.,class=multipole;
 esave,  file="temp/MB.errors";
 system, "slhc/errors/corr_MB_phaseI";
 call,   file="temp/MB_corr_setting.mad";
!exec reset_MB_corr;

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
! correction of triplet and D1 (only in collision, not compatible SLHC)
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if (NRJ>5000.0000)
{
  option, echo, warn, -info;
 select, flag=error, clear;
 select, flag=error, pattern=MQX,  class=multipole;
 select, flag=error, pattern=MBX,  class=multipole;
 select, flag=error, pattern=MBRC, class=multipole;
 esave,  file="temp/tripD1D2.errors";
 system, "db5/toolkit/corr_tripD1";
 call,   file="temp/MCX_setting.mad";
};


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!    correct orbit distortion resulting from other magnets
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if((ON_A1S)^2+(ON_A1r)^2+(ON_B1S)^2+(ON_B1r)^2 >0)
{
 exec, initial_micado(4);
 exec, initial_micado(4);
 exec, initial_micado(4);
 exec, initial_micado(4);
 };
 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!    final orbit correction before applying crossing scheme
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if((ON_A1S)^2+(ON_A1r)^2+(ON_B1S)^2+(ON_B1r)^2 >0)
{ exec, final_micado(0.004); };


! end of orbit correction, now switch ON crossing scheme
 on_x1:=1; on_sep1:=1;
 on_x2:=1; on_sep2:=1; on_alice:=1;
 on_x5:=1; on_sep5:=1;
 on_x8:=1; on_sep8:=1; on_lhcb :=-1;
 
! coguess,x=x.ip1,px=px.ip1,y=y.ip1,py=py.ip1;

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!   Fine tuning of coupling after CO correction and with Xscheme
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

call,file="slhc/errors/FineCouplingCorrection.madx";

 
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                 matching of tune and chromaticity
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
if(mylhcbeam==1)
{kqtf=kqtf.b1;kqtd=kqtd.b1;kqtf.b1:=kqtf;kqtd.b1:=kqtd;
 ksf=ksf.b1;ksd=ksd.b1;ksf.b1:=ksf;ksd.b1:=ksd;};

if(mylhcbeam>1)
{kqtf=kqtf.b2;kqtd=kqtd.b2;kqtf.b2:=kqtf;kqtd.b2:=kqtd;
 ksf=ksf.b2;ksd=ksd.b2;ksf.b2:=ksf;ksd.b2:=ksd;};

 match;
 global, q1=qx0, q2=qy0;
 vary,   name=kqtf, step=1.0E-7 ;
 vary,   name=kqtd, step=1.0E-7 ;
 lmdif,  calls=100, tolerance=1.0E-21;
 endmatch;
 
 match,chrom;
 global, dq1=qprime, dq2=qprime;
 vary,   name=ksf;
 vary,   name=ksd;
 lmdif,  calls=100, tolerance=1.0E-21;
 endmatch; 
  
 match,chrom;
 global, dq1=qprime, dq2=qprime;
 global, q1=qx0, q2=qy0;
 vary,   name=ksf;
 vary,   name=ksd;
 vary,   name=kqtf, step=1.0E-7 ;
 vary,   name=kqtd, step=1.0E-7 ;
 lmdif,  calls=500, tolerance=1.0E-21;
 endmatch;


!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                 Check corrector strength
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
nrj=NRJ*1.E9;
 if (mylhcbeam == 1){call, file="db5/toolkit/corr_value_b1.madx";}
 if (mylhcbeam >  1){call, file="db5/toolkit/corr_value_b2.madx";}

!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!                 final twiss before sending to sixtrack
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

/*
twiss;
 value, table(summ,q1);
 value, table(summ,q2);
 value, table(summ,dq1);
 value, table(summ,dq2);
 sixtrack, radius=0.017;
*/

!/*
if (NRJ<4999.9999) {VRF400:=8. ;LAGRF400.B1=0.5;LAGRF400.B2=0.;};
if (NRJ>5000.0000) {VRF400:=16.;LAGRF400.B1=0.5;LAGRF400.B2=0.;};
twiss;
!emit;
 value, table(summ,q1);
 value, table(summ,q2);
 value, table(summ,dq1);
 value, table(summ,dq2);
 sixtrack,cavall, radius=0.017;
!*/ 

select, flag=twiss, clear;
twiss,file="last_twiss.%SEEDRAN";
System,"gzip last_twiss.%SEEDRAN";

stop;
