CorrectionChain

0001 function [...
0002     rcor,...            % corrected lattice
0003     ch,...              % final H cor values
0004     cv,...              % final V cor values
0005     cq,...              % final Quad cor values
0006     cs,...              % final Skew Quad cor values
0007     inCOD,...
0008     d0,de,dc...         % lattice data structures d0=no err, de=err, dc=cor
0009     ]=CorrectionChain(...
0010     rerr,...            %1  initial lattice
0011     r0,...              %2  model lattice
0012     indBPM,...          %3  bpm index
0013     indHCor,...         %4  h steerers index
0014     indVCor,...         %5  v steerers index
0015     indSkewQuadCor,...  %6  skew quad index
0016     indQuadCor,...      %7  quadrupole correctors index
0017     Neig,...            %8  number of eigen vectors [NeigorbitH, NeigorbitV, NeigQuadrdt, Neigdispv, Neigdisph,neig rdt corr, SkewQuadRDT]
0018     corrorder,...       %9  correction order 1: orbit, 2: tune, 3: skewquad disp v 4: quad disp h 5: quad RDT 6: skew RDT
0019     ModelRM,...          %10 response matrices
0020     speclab,...          %11 response matrices
0021     verbose)             %12 verbose (false): if true print out all relevat quantities after each step in corrorder
0022 %
0023 %
0024 % performs a loop of corrections as described in corparam
0025 %
0026 % [...
0027 %     rcor,...            % corrected lattice
0028 %     ch,...              % final H cor values
0029 %     cv,...              % final V cor values
0030 %     cq,...              % final Quad cor values
0031 %     cs...               % final Skew Quad cor valu
0032 %     ]=CorrectionChain(...
0033 %     rerr,...            %1  initial lattice
0034 %     r0,...              %2  model lattice
0035 %     indBPM,...          %3  bpm index
0036 %     indHCor,...         %4  h steerers index
0037 %     indVCor,...         %5  v steerers index
0038 %     indSkewQuadCor,...  %6  skew quad index
0039 %     indQuadCor,...      %7  quadrupole correctors index
0040 %     Neig,...            %8  number of eigen vectors [NeigorbitH, NeigorbitV, NeigQuadrdt, Neigdispv, Neigdisph,neig rdt corr, SkewQuadRDT]
0041 %     corrorder,...       %9  correction order 1: orbit, 2: tune, 3: skewquad disp v 4: quad disp h 5: quad RDT 6: skew RDT
0042 %     ModelRM,...          %10 response matrices
0043 %     speclab,...          %11 response matrices
0044 %     verbose)             %12 verbose (false): if true print out all relevat quantities after each step in corrorder
0045 %
0046 %
0047 %
0048 % number of eigenvalues for SVD
0049 % neigSteererH     = Neig(1); % H orbit correction
0050 % neigSteererV     = Neig(2); % V orbit correction
0051 % neigSkew         = Neig(3); % skew quadrupoles vert. disp. correction
0052 % neigQuad         = Neig(4); % quadrupoles hor. disp. correction, beta correction and phase correction
0053 % neigQuadFit      = Neig(5); % quadrupole errors fit
0054 % neigDipFit       = Neig(6); % dipole errors fit
0055 % neigSkewFit      = Neig(7); % skew quad fit
0056 %
0057 %
0058 % corrorder=[0:7];
0059 %
0060 %     '(-1 ): RF cavity frequency and time lag tuning '...
0061 %     '( 0): open trajectory (finds closed orbit) '...
0062 %     '( 1): orbit '...
0063 %     '( 2): tune '...
0064 %     '( 3): chromaticity '...
0065 %     '( 4): dispersion '...
0066 %     '( 5): dispersion free steering '...
0067 %     '( 6): rdt + dispersion correction '...
0068 %     '( 7): fit errors model and correct model quad RDT + dispersion (6) '
0069 %
0070 % if the correction fails, nothing is done to the lattice.
0071 %
0072 %see also:
0073 % findrespmat
0074 % qemsvd_mod
0075 % atsetRFCavity
0076 % atfirstturntrajectory
0077 % atcorrectorbit
0078 % fittunedelta2fam
0079 % atmatchchromdelta
0080 % atcorrectdispersion
0081 % atdispersionfreesteering
0082 % atRDTdispersioncorrection
0083 % FitResponseMatrixAndDispersion
0084 %
0085 
0086 
0087 t0=tic;
0088 disp('>>>>>  Correction sequence started  <<<<<')
0089 Nbpm=length(indBPM);
0090 NVcor=length(indVCor);
0091 NHcor=length(indHCor);
0092 NScor=length(indSkewQuadCor);
0093 NQcor=length(indQuadCor);
0094 disp(['     # BPM  : ' num2str(Nbpm) ]);
0095 disp(['     # V cor: ' num2str(NVcor)]);
0096 disp(['     # H cor: ' num2str(NHcor)]);
0097 disp(['     # S cor: ' num2str(NScor)]);
0098 disp(['     # Q cor: ' num2str(NQcor)]);
0099 disp('>>>>>  --------------------------  <<<<<')
0100 
0101 if nargin<12
0102     verbose=false;
0103 end
0104 
0105 if nargin<11
0106     speclab='';
0107 end
0108 
0109 if nargin<10
0110     ModelRM=[];
0111 else
0112     disp(' --- - - - - - - - - - - - - - - - - - - ---')
0113     disp(' ---                                     ---')
0114     disp(' --- FAST MODE: MODEL RM ARE BEING USED! ---')
0115     disp(' ---                                     ---')
0116     disp(' --- - - - - - - - - - - - - - - - - - - ---')
0117 end
0118 
0119 
0120 if nargin<9
0121     
0122     corrorder=[0:3];
0123     % 0: open trajectory
0124     % 1: orbit
0125     % 2: tune
0126     % 3: chromaticity
0127     
0128 end
0129 
0130 if nargin<8  % default number of eigenvectors is 100
0131     disp('100 eigen vectors for all corrections')
0132     Neig = 100 * ones(7,1);
0133 end
0134 % number of eigenvalues for SVD
0135 neigSteererH     =Neig(1); % H orbit correction
0136 neigSteererV     =Neig(2); % V orbit correction
0137 neigSkew         =Neig(3); % skew quadrupoles vert. disp. correction
0138 neigQuad         =Neig(4); % quadrupoles hor. disp. correction, beta correction and phase correction
0139 neigQuadFit      =Neig(5); % quadrupole errors fit
0140 neigDipFit       =Neig(6); % dipole errors fit
0141 neigSkewFit      =Neig(7); % skew quad fit
0142 
0143 ch=atgetfieldvalues(rerr,indHCor,'PolynomB',{1,1});
0144 cv=atgetfieldvalues(rerr,indVCor,'PolynomA',{1,1});
0145 cq=atgetfieldvalues(rerr,indQuadCor,'PolynomB',{1,2});
0146 cs=atgetfieldvalues(rerr,indSkewQuadCor,'PolynomA',{1,2});
0147 
0148 rerrINIT=rerr; % initial lattice with errors to compute correction (PolynomB(2) stores everything!)
0149 
0150 % display selected correction order
0151 disp('>>>>>      Correction STRATEGY :       <<<<<')
0152 disp('.')
0153 for iorddisp=corrorder
0154     switch iorddisp
0155         case -1
0156             disp('          RF cavity')
0157         case 0
0158             disp('          open trajectory (steerers)')
0159         case 1
0160             disp('          orbit (steerers) ')
0161         case 2
0162             disp('          tune (quadrupoles, 2 families)')
0163         case 3
0164             disp('          chromaticity (sextupoles, 2 families)')
0165         case 4
0166             disp('          dispersion (quadrupoles)')
0167         case 5
0168             disp('          dispersion free steering (correctors)')
0169         case 6
0170             disp('          RDT + dispersion (quadrupoles) ')
0171         case 7
0172             disp('          Fit Quad+Dip Errors ')
0173             disp('          Correct RDT and Dispersion of fitted model ')
0174         otherwise
0175     end
0176 end
0177 disp('.')
0178 disp('>>>>>  --------------------------  <<<<<')
0179 
0180 
0181 inCOD=[0 0 0 0 0 0]';
0182 [l,~,~]=atlinopt(r0,0,indBPM);
0183 refdispersion=zeros(2,length(indBPM));
0184 refdispersion(1,:)=arrayfun(@(a)a.Dispersion(1),l);
0185 refdispersion(2,:)=arrayfun(@(a)a.Dispersion(3),l);
0186 
0187 % perform correction
0188 iicor=1;
0189 for cor=corrorder
0190     tic;
0191     mesgcor=['Correction Step: ' num2str(iicor) '/' num2str(length(corrorder))];
0192     disp(mesgcor);
0193     iicor=iicor+1;
0194     
0195     rerr0=rerr;% for correction display.
0196     inCODe=inCOD;
0197     
0198     switch cor
0199         
0200         case -1
0201             %% set rfcavity
0202             
0203             % decide if radiation is on or off in the lattice by looking
0204             % for Rad in Pass Methods
0205             radon=any(cellfun(@(a)sum(ismember('Rad',a.PassMethod))==3,rerr));
0206             
0207             % get cavity settings
0208             indrfc=find(atgetcells(rerr,'Frequency'));
0209             rfv=sum(atgetfieldvalues(rerr,indrfc,'Voltage'));
0210             harm=atgetfieldvalues(rerr,indrfc(1),'HarmNumber');
0211             tlag=atgetfieldvalues(rerr,indrfc(1),'TimeLag');
0212             
0213             disp(['Set RF cavity. '....
0214                 num2str(rfv*1e-6) ' MV, '...
0215                 num2str(harm) ' buckets, '...
0216                 num2str(radon) ' radiation '...
0217                 ]);
0218             
0219             % set cavities rf frequency and time lag for lattice with
0220             % errors
0221           
0222             rerr=atsetRFCavityErr(rerr,rfv,radon,harm,inCOD);
0223             
0224             [...
0225                 rerr,...
0226                 inCOD...
0227                 ]=atRFcorrection(...
0228                 rerr,...
0229                 indBPM,...
0230                 inCOD,...
0231                 [1 1 1 1 1],...
0232                 1,...
0233                 ModelRM);
0234 
0235         case 0
0236             %% OPEN TRAJECTORY
0237             excursion=1e-3;
0238             disp(['Open Trajectory Correction started. '....
0239                 num2str(excursion*1e3) 'mm escursion accepted']);
0240             
0241             if isempty(ModelRM)
0242                 ModelRMtr=r0;
0243             else
0244                 ModelRMtr=ModelRM;
0245             end
0246             [rerr,inCOD]=atfirstturntrajectory(...
0247                 rerr,...
0248                 inCOD,...
0249                 indBPM,...
0250                 indHCor,...
0251                 indVCor,...
0252                 excursion,...
0253                 200,...
0254                 [false true],...
0255                 ModelRMtr);
0256             
0257             ch=atgetfieldvalues(rerr,indHCor,'PolynomB',{1,1});
0258             cv=atgetfieldvalues(rerr,indVCor,'PolynomA',{1,1});
0259             
0260             try
0261             catch exc
0262                 
0263                 getReport(exc,'extended');
0264                 error('Failed trajectory correction');
0265                 
0266             end
0267             
0268         case 1
0269             %% ORBIT CORRECTION
0270             disp(['Steerers to fix orbit H ' ...
0271                 ' using ' num2str(neigSteererH) ' eig']);
0272             disp(['Steerers to fix orbit V ' ...
0273                 ' using ' num2str(neigSteererV) ' eig']);
0274             
0275             
0276             [rerr,inCOD]=atcorrectorbit(rerr,...
0277                 indBPM,...
0278                 indHCor,...
0279                 indVCor,...
0280                 inCOD,...
0281                 [[floor(linspace(1,neigSteererH,10)),neigSteererH,neigSteererH];...
0282                 [floor(linspace(1,neigSteererV,10)),neigSteererV,neigSteererV]]',...
0283                 [false true],... dpp correction and mean to zero
0284                 1.0,...
0285                 ModelRM);
0286                         
0287             
0288             ch=atgetfieldvalues(rerr,indHCor,'PolynomB',{1,1});
0289             cv=atgetfieldvalues(rerr,indVCor,'PolynomA',{1,1});
0290             
0291             try   
0292             catch exc
0293                 
0294                 getReport(exc);
0295                 warning('Failed Orbit correction')
0296                 
0297             end
0298         case 2
0299             %% TUNE MATCH (2 families)
0300             
0301             try % tune rematch
0302                 
0303                 rerr=fittunedelta2fam(rerr,r0);
0304                 
0305             catch exc
0306                 
0307                 getReport(exc);
0308                 save('latticeFailedTunecor.mat','rerr','r0');
0309                 disp('Could not match Tune');
0310                 warning('Could not match Tune');
0311                 
0312             end
0313             
0314         case 3
0315             %% chromaticity correction
0316             
0317             disp(' - - - -  chromaticty correction - - - - - ');
0318             disp('All SF and All SD moved by a constant value');
0319             
0320             try
0321                 
0322                 indS=find(atgetcells(r0,'Class','Sextupole'))';
0323                 pbsxt=atgetfieldvalues(r0,indS,'PolynomB',{1,3});
0324                 indSF=indS(pbsxt>0);
0325                 indSD=indS(pbsxt<0);
0326                 
0327                 [~,~,chrom]=atlinopt(r0,0,1);disp(['Nominal chrom: ' num2str(chrom,'%2.3f, ')]);
0328                 [~,~,chrome]=atlinopt(rerr,0,1);disp(['Initial chrom: ' num2str(chrome,'%2.3f, ')]);
0329                 
0330                 rerr=atmatchchromdelta(rerr,chrom,{indSF,indSD});
0331                 
0332                 [~,~,chromcor]=atlinopt(rerr,0,1);disp(['Final chrom: ' num2str(chromcor,'%2.3f, ')]);
0333                 
0334                 
0335             catch exc
0336                 
0337                 getReport(exc);
0338                 warning('Failed chromaticty correction')
0339                 
0340             end
0341             
0342         case 4
0343             %% DISPERSION CORRECTION
0344             disp(['Normal quadrupoles to fix dispersion H ' ...
0345                 ' using ' num2str(neigQuad) ' eig']);
0346             disp(['Skew quadrupoles to fix dispersion V ' ...
0347                 ' using ' num2str(neigSkew) ' eig']);
0348             
0349             
0350             [rerr,inCOD]=atcorrectdispersion(rerr,...
0351                 indBPM,...
0352                 indQuadCor,...
0353                 indSkewQuadCor,...
0354                 inCOD,...
0355                 [[floor(linspace(20,neigQuad,5)),neigQuad,neigQuad];...
0356                 [floor(linspace(20,neigSkew,5)),neigSkew,neigSkew]]',...
0357                 [true true],... dpp correction and mean to zero
0358                 1.0,...
0359                 ModelRM,...
0360                 refdispersion,...
0361                 [],...
0362                 true);
0363             
0364             
0365             cq=atgetfieldvalues(rerr,indQuadCor,'PolynomB',{1,2});
0366             cs=atgetfieldvalues(rerr,indSkewQuadCor,'PolynomA',{1,2});
0367             
0368             try
0369                 
0370             catch exc
0371                 
0372                 getReport(exc);
0373                 warning('Failed Orbit correction')
0374                 
0375             end
0376             
0377         case 5
0378             %% ORBIT+DISPERSION CORRECTION (dispersion free steering)
0379             disp(['Steerers to fix orbit and dispersion H ' ...
0380                 ' using ' num2str(neigSteererH) ' eig']);
0381             disp(['Steerers to fix orbit and dispersion V ' ...
0382                 ' using ' num2str(neigSteererV) ' eig']);
0383             
0384             
0385             [rerr,inCOD]=atdispersionfreesteering(...
0386                 rerr,...
0387                 indBPM,...
0388                 indHCor,...
0389                 indVCor,...
0390                 inCOD,...
0391                 [[floor(linspace(20,neigSteererH,10)),neigSteererH,neigSteererH];...
0392                 [floor(linspace(20,neigSteererV,10)),neigSteererV,neigSteererV]]',...
0393                 [true true],...
0394                 1.0,...
0395                 0.9,...
0396                 ModelRM,...
0397                 zeros(2,length(indBPM)),...
0398                 refdispersion,...
0399                 [],...
0400                 true);
0401             
0402             ch=atgetfieldvalues(rerr,indHCor,'PolynomB',{1,1});
0403             cv=atgetfieldvalues(rerr,indVCor,'PolynomA',{1,1});
0404             
0405             try
0406                 
0407             catch exc
0408                 
0409                 getReport(exc);
0410                 warning('Failed Orbit-dispersion (DFS) correction')
0411                 
0412             end
0413         case 6
0414             %% RDT+DISPERSION CORRECTION
0415             disp(['Quadrupoles to fix RDT, tune and dispersion H ' ...
0416                 ' using ' num2str(neigQuad) ' eig']);
0417             disp(['Steerers to fix RDT and dispersion V ' ...
0418                 ' using ' num2str(neigSkew) ' eig']);
0419             
0420             [rerr,inCOD]=atRDTdispersioncorrection(...
0421                 rerr,...
0422                 r0,...
0423                 indBPM,...
0424                 indQuadCor,...
0425                 indSkewQuadCor,...
0426                 inCOD,...
0427                 [[floor(linspace(1,neigQuad,5)),neigQuad,neigQuad];...
0428                 [floor(linspace(1,neigSkew,5)),neigSkew,neigSkew]]',...
0429                 [false],...
0430                 1.0,...
0431                 [0.8 0.1 0.8],...
0432                 ModelRM);
0433             
0434             cq=atgetfieldvalues(rerr,indQuadCor,'PolynomB',{1,2});
0435             cs=atgetfieldvalues(rerr,indSkewQuadCor,'PolynomA',{1,2});
0436             
0437             try
0438                 
0439             catch exc
0440                 
0441                 getReport(exc);
0442                 warning('Failed RDT and dispersion correction')
0443                 
0444             end
0445             
0446        
0447         case 7
0448             %% RDT+DISPERSION CORRECTION from lattice error model
0449                 % fit lattice errors model
0450                 [rfit]=FitResponseMatrixAndDispersion(...
0451                     rerr,...
0452                     r0,...
0453                     inCOD,...
0454                     indBPM,...
0455                     indHCor(1:9*2:end),... % 4 correctors, 1 every 8 cells
0456                     indHCor(1:9*2:end),...  % 4 correctors, 1 every 8 cells
0457                     [neigQuadFit,neigDipFit,neigSkewFit,neigDipFit],...
0458                     4,...
0459                     [speclab 'fitrm']);
0460                 
0461                 % get change of strength of correctors
0462                 fq=atgetfieldvalues(rfit,indQuadCor,'PolynomB',{1,2});
0463                 fs=atgetfieldvalues(rfit,indSkewQuadCor,'PolynomA',{1,2});
0464                 
0465                 % correct RDT and dispersion of fitted error model
0466                 [~,inCOD,fcq,fcs]=atRDTdispersioncorrection(...
0467                     rfit,... <<--- fitted error model! not lattice with errors!
0468                     r0,...
0469                     indBPM,...
0470                     indQuadCor,...
0471                     indSkewQuadCor,...
0472                     inCOD,...
0473                     [[floor(linspace(1,neigQuad,5)),neigQuad,neigQuad];...
0474                     [floor(linspace(1,neigSkew,5)),neigSkew,neigSkew]]',...
0475                     [true],...
0476                     1.0,...
0477                     [0.8 0.1 0.8],...
0478                     ModelRM);
0479                 
0480                 %fcq=atgetfieldvalues(rfitcor,indQuadCor,'PolynomB',{1,2});
0481                 %fcs=atgetfieldvalues(rfitcor,indSkewQuadCor,'PolynomA',{1,2});
0482                 
0483                 % store proposed correction
0484                 dcq(1,:)=(fcq-fq);
0485                 dcs(1,:)=(fcs-fs);
0486                 
0487             
0488             % set delta correctors strength in lattice with errors.
0489             cq=atgetfieldvalues(rerr,indQuadCor,'PolynomB',{1,2});
0490             cs=atgetfieldvalues(rerr,indSkewQuadCor,'PolynomA',{1,2});
0491             
0492             cq=cq+dcq'; %  add proposed correction on fitted lattice
0493             cs=cs+dcs';
0494             rerr=atsetfieldvalues(rerr,indQuadCor,'PolynomB',{1,2},cq);
0495             rerr=atsetfieldvalues(rerr,indSkewQuadCor,'PolynomA',{1,2},cs);
0496             
0497             
0498             try     
0499             catch exc
0500                 
0501                 getReport(exc);
0502                 warning('Failed model fit and RDT + dispersion correction')
0503                 
0504             end
0505             
0506             
0507         
0508         otherwise
0509             warning([num2str(cor) ': not a possible correction. [0:7]: '...
0510                 '( 0): open trajectory (finds closed orbit) '...
0511                 '( 1): orbit '...
0512                 '( 2): tune '...
0513                 '( 3): chromaticity '...
0514                 '( 4): dispersion '...
0515                 '( 5): dispersion free steering '...
0516                 '( 6): rdt + dispersion correction '...
0517                 '( 7): fit errors model and correct model quad RDT + dispersion (6) '])
0518     end
0519     
0520     %     %% set corrector multipoles
0521     %     rerr=SetCorMult(rerr);
0522     %
0523     
0524     % %%apply correctors PS limits
0525     %rerr=SetCorLimits(rerr);
0526     if verbose
0527         %% display correction effect at every step
0528         DisplayCorrectionEffect(r0,rerr0,rerr,inCODe,inCOD,1:length(r0),indHCor,indVCor,indQuadCor,indSkewQuadCor);
0529     end
0530     
0531     disp(['Finished: ' mesgcor])
0532     toc;
0533 end
0534 
0535 d0=[];de=[];dc=[];
0536 if ~verbose
0537     %% display correction effect from begining to end.
0538  
0539     [d0,de,dc]=DisplayCorrectionEffect(r0,rerrINIT,rerr,inCODe,inCOD,[1:length(r0)]',indHCor,indVCor,indQuadCor,indSkewQuadCor);
0540 end
0541 rcor=rerr;
0542 
0543 tend=toc(t0);
0544 disp(['Time for correction chain: ' num2str((tend-t0)/60) ' minutes'])
0545 return