From aa371423c7b38c84f4881101d9adec23d20e7d97 Mon Sep 17 00:00:00 2001
From: Laurent Farvacque <laurent.farvacque@gmail.com>
Date: Fri, 27 Dec 2019 16:51:01 +0100
Subject: [PATCH 1/4] added DQs and SFs to the error fit

---
 qem/qempanelfitq.m         |  9 ++++--
 qem/qempanelfitqAnalytic.m | 10 +++---
 qem/qempanelset.m          | 66 ++++++++++----------------------------
 3 files changed, 29 insertions(+), 56 deletions(-)

diff --git a/qem/qempanelfitq.m b/qem/qempanelfitq.m
index 8deddd3d..0e2728ad 100644
--- a/qem/qempanelfitq.m
+++ b/qem/qempanelfitq.m
@@ -16,14 +16,16 @@ function newkn=qempanelfitq(qemb,qemres,semres,okfit,okbpm,varargin)
 %
 %KN:        New quadrupole strengths
 
+global GDRESP
+
     function elem=setk(elem,dk)
         strength=elem.PolynomB(2)+dk;
         elem.K=strength;
         elem.PolynomB(2)=strength;
     end
 
-if nargin<5, okbpm=[]; end;
-if nargin<4, okfit=[]; end;
+if nargin<5, okbpm=[]; end
+if nargin<4, okfit=[]; end
 narg=length(varargin);
 if narg>0 && isa(varargin{narg},'function_handle')
     disparg={@(i,itot) nselect(2,varargin{narg},i,itot)};
@@ -33,7 +35,7 @@ else
 end
 if narg>0 && isstruct(varargin{narg})
     mode=varargin{narg};
-    narg=narg-1;
+    narg=narg-1; %#ok<NASGU>
 else
     mode=struct;
 end
@@ -47,6 +49,7 @@ nvst=length(svidx);
 dpp=o0(5);
 [~,fractunes]=atlinopt(qemb.at,dpp,[],o0);
 dresp=qemderiv(qemb.at,dpp,@setk,0.0001,qemres.qpidx(okfit),qemres.bpmidx,shidx,svidx,disparg{:});
+GDRESP=dresp;
 
 [rh,rh2v,rv2h,rv,frh,frv]=qemdecode(qemb.at,qemres.ct,qemres,...
     qemres.resph,semres.respv,semres.resph,qemres.respv,...
diff --git a/qem/qempanelfitqAnalytic.m b/qem/qempanelfitqAnalytic.m
index 4ea0c77c..b6043b86 100644
--- a/qem/qempanelfitqAnalytic.m
+++ b/qem/qempanelfitqAnalytic.m
@@ -16,6 +16,8 @@ function newkn=qempanelfitqAnalytic(qemb,qemres,semres,okfit,okbpm,varargin)
 %
 %KN:        New quadrupole strengths
 
+global GDRESPA
+
 disp('ANALYTIC quadrupole fit QEMPANEL')
 
     function elem=setk(elem,dk)
@@ -24,8 +26,8 @@ disp('ANALYTIC quadrupole fit QEMPANEL')
         elem.PolynomB(2)=strength;
     end
 
-if nargin<5, okbpm=[]; end;
-if nargin<4, okfit=[]; end;
+if nargin<5, okbpm=[]; end
+if nargin<4, okfit=[]; end
 narg=length(varargin);
 if narg>0 && isa(varargin{narg},'function_handle')
     disparg={@(i,itot) nselect(2,varargin{narg},i,itot)};
@@ -35,7 +37,7 @@ else
 end
 if narg>0 && isstruct(varargin{narg})
     mode=varargin{narg};
-    narg=narg-1;
+    narg=narg-1; %#ok<NASGU>
 else
     mode=struct;
 end
@@ -55,7 +57,7 @@ dpp=o0(5);
 %dresp=qemderivAnalyticDebuggingVersion(qemb.at,dpp,@setk,0.0001,qemres.qpidx(okfit),qemres.bpmidx,shidx,svidx,disparg{:});
 dresp=qemderivAnalytic(qemb.at,dpp,@setk,0.0001,qemres.qpidx(okfit),qemres.bpmidx,shidx,svidx,disparg{:});
 %drespN=qemderiv(qemb.at,dpp,@setk,0.0001,qemres.qpidx(okfit),qemres.bpmidx,shidx,svidx,disparg{:});
-
+GDRESPA=dresp;
 % drespA = dresp;
 
 % % % lines for tests 
diff --git a/qem/qempanelset.m b/qem/qempanelset.m
index 619cfe40..a67001f9 100644
--- a/qem/qempanelset.m
+++ b/qem/qempanelset.m
@@ -35,7 +35,6 @@ qemres.opticsdir=fullfile(dirname,'optics','');
 [~,theodir]=system(['readlink ' fullfile(dirname,'optics','')]);
 
 try
-    
     simu = tango.Device('sys/ringsimulator/ebs');
     qemres.simufile=simu.get_property('RingFile');
     
@@ -55,11 +54,10 @@ warning('check SB are here')
 atmod=sr.model(ring,'reduce',true); % for number of DQ (96, not 131 splitted) for example
 atmodpinholes=sr.model(sr.model.addpinholes(atmod.ring)); % this lattice has pinhole cameras but DL & DQ indexes are bad
 
-indhcor = sort([atmod.get(0,'steerhdq')]);
-indvcor = sort([atmod.get(0,'steerv')]);
+indhcor = sort(atmod.get(0,'steerhdq'));
+indvcor = sort(atmod.get(0,'steerv'));
 
 if full_partial
-    
      qemres.hlist=1:12:length(indhcor);
      qemres.hlist(1) = 4;
      qemres.vlist=1:9:length(indvcor);
@@ -133,33 +131,6 @@ atmod.settune(qemres.fulltunes,'QF1QD2');
 
 mach = atmod.ring;
 
-try
-    % % note: quadrupole and skew quaadrupole strengths are assigned by qemat.m
-    
-    % set sextupole correction strengths from accelerator to AT lattice
-    indsext=atmod.get(0,'sx')'; %get sextupole indexes
-    
-    Ksext=sr.read_mag(sr.sxname([1:length(indsext)]));%,'Tango://ebs-simu:10000/srmag/m-s')');
-    disp(['' num2str(length(isnan(Ksext))) ...
-        ' sextupoles corrections std: ' num2str(std(Ksext))]);
-    
-    Ksext0 = atgetfieldvalues(mach,indsext(~isnan(Ksext)),'PolynomB',{1,3});
-    
-    mach=atsetfieldvalues(mach,indsext(~isnan(Ksext)),'PolynomB',{1,3},Ksext0(~isnan(Ksext))+Ksext(~isnan(Ksext))');
-     
-    % set octupole correction strengths from accelerator to AT lattice
-    indoct=atmod.get(0,'oc')'; %get sextupole indexes
-    Koct=sr.read_mag(sr.ocname([1:length(indoct)]));%,'Tango://ebs-simu:10000/srmag/m-s')');
-    disp(['' num2str(length(isnan(Koct))) ...
-        ' octupole corrections std: ' num2str(std(Koct))]);
-    Koct0 = atgetfieldvalues(mach,indoct(~isnan(Ksext)),'PolynomB',{1,4});
-    
-    mach=atsetfieldvalues(mach,indoct(~isnan(Koct)),'PolynomB',{1,4},Koct0(~isnan(Ksext))+Koct(~isnan(Koct))');
-    
-catch
-    warning('qem:noinit','Could not read sextupoles values from S.R. . KEPT as in reference optics, no corrections.');
-end
-
 % store lattice in qemres model
 qemres.at=mach;
 
@@ -182,34 +153,32 @@ qemres.fractunes=periods*tunes-fix(periods*tunes);
 
 % this first list is ok for dipoles, not ok for skew quadrupoles!
 % qemres.dipidx=find(atgetcells(atmod.ring,'FamName','[DJ]L[12][AE]_[15]','DQ1[B]\w*','DQ2C\w*'))';%,'DQ1[BD]' atmod.get(0,'di');
-dlind=find(atgetcells(atmod.ring,'FamName','[DJ]L[12][ABDE]_[15]\w*'))';%,'DQ1[BD]' atmod.get(0,'di');
-mlind=find(atgetcells(atmod.ring,'FamName','ml\w*'))';%,'DQ1[BD]' atmod.get(0,'di');
-qemres.dipidx=sort([dlind,mlind(1:2:end)]);%,'DQ1[BD]' atmod.get(0,'di');
-bang=atmod.getfieldvalue(qemres.dipidx,'BendingAngle');
-qemres.dipidx=qemres.dipidx(bang>0.00001);
-%qemres.dipidx=qemres.dipidx(1);
-
+dlmask=atgetcells(atmod.ring,'FamName','[DJ]L[12][ABDE]_5\w*');
+dqmask=atmod.select('dq');
+qemres.dipidx=find(dlmask | dqmask);
 
 % AnalyticResponseMatrix.m must be modified to include DQ! also qemderiv*.m
 % line 96 should be modified to get this working
 
- qemres.qpidx= sort([atmod.get(0,'qp');])'; % atmod.get(0,'dq') fit locations
-% qemres.qpidx= [atmod.get(0,'qp');]'; % fit locations
+qpmask=atmod.select('qp');
+sfmask=atgetcells(atmod.ring,'FamName','S[FIJ]2[AE]');
+qf4mask=atgetcells(atmod.ring,'FamName','QF4[ABDEIJ]');
+fmask=qpmask | dqmask | sfmask;  % All QPs, DQs, SF2s
+
+qemres.qpidx=find(fmask)';  % All QPs, DQs, SF2s
 
 % % warning('only quadrupoles, no DQ!')
-qemres.quadidx = atmod.get(0,'qp')' ; % only quadrupoles
+qemres.quadidx = find(dqmask);	% Only quadrupoles
 
-% qemres.dqidx=atmod.get(0,'dq')';
-qemres.dqidx=[atmod.get(0,'dq');]';
+qemres.dqidx=find(dqmask);
 
-qemres.qcoridx=atmod.get(0,'qp')'; % either quadrupoles or DQ 
+qemres.qcoridx=find(qpmask)';    % Only quadrupoles
 qemres.qcorl=atgetfieldvalues(qemres.at(qemres.qcoridx),'Length');
 
-%qemres.qpfit = ismember(qemres.qpidx,qemres.qcoridx);%keep all quads no DQ
-qemres.qpfit = true(size(qemres.qpidx));  % keep all quads and DQ
+qemres.qpfit = fmask(fmask) & ~qf4mask(fmask);
 qemres.ql=atgetfieldvalues(qemres.at(qemres.qpidx),'Length');
 
-qemres.sextidx=atmod.get(0,'sx')';
+qemres.sextidx=find(atmod.select('sx'));
 sxl=atgetfieldvalues(qemres.at(qemres.sextidx),'Length');
 sxl(~(sxl~=0))=1;   % Handles NaN correctly
 qemres.sxl=sxl;
@@ -316,8 +285,7 @@ if nargout >= 3
     end
     
     Nq=textscan(fopen(fullfile(qemres.datadir,'quadcor.dat'),'r'),'%s %f','HeaderLines',7);
-    allqdq=sort(unique([qemres.quadidx,qemres.dqidx([true,diff(qemres.dqidx)>2])]));
-    [~,indq]=ismember(qemres.qcoridx,allqdq);% quad in quad dq list
+    indq=qpmask(qpmask | dqmask);
     corq = Nq{2}(indq); %select Quadrupoles
     qemb(2).cor = corq;
     Sq=textscan(fopen(fullfile(qemres.datadir,'skewcor.dat'),'r'),'%s %f','HeaderLines',7);
-- 
GitLab


From 23bffe619d186065e21911cf44c04a979b1d8300 Mon Sep 17 00:00:00 2001
From: Laurent Farvacque <laurent.farvacque@esrf.fr>
Date: Thu, 2 Jan 2020 16:49:45 +0100
Subject: [PATCH 2/4] absolute/relative

---
 +sr/fold.m         | 37 +++++++++++++++---------
 qem/qcheck.m       | 48 +++++++++++++++++--------------
 qem/qempaneldisp.m | 72 +++++++++++++++++++++++++++-------------------
 qem/qempanelset.m  | 52 +++++++++++++++++++--------------
 qem/qemresp.m      |  1 -
 sem/semtilt.m      |  2 +-
 6 files changed, 124 insertions(+), 88 deletions(-)

diff --git a/+sr/fold.m b/+sr/fold.m
index 63c2744f..196fcefa 100644
--- a/+sr/fold.m
+++ b/+sr/fold.m
@@ -1,21 +1,30 @@
 function r=fold(a,dim)
 %A=SRFOLD(B)			reshapes data for SR symmetry
 %
-%	A(32*n,1) -> R{n,32}
-
-
-r{1}=[a(1),a(end)];
-r{2}=[a(2),a(end-1)];
-
-a=a(3:end-2);
-
-nq=numel(a)/32;
-
-for iq=1:nq
-    r{2+iq}=a(iq:nq:end);
+%	A(64*n,1) -> R{n,64}
+
+if isvector(a)
+    if numel(a) == 514      % Special case of quadrupoles
+        r=[sr.fold(a(2:end-1));NaN*ones(3,64)];
+        r(1:2,1)=NaN;
+        r(1:2,64)=NaN;
+        r(9:11,1)=a(1:3);
+        r(11:-1:9,64)=a(end-2:end);
+    else
+        n2=numel(a)/32;
+        b=reshape(a(:),n2,32);
+        if bitand(n2,1) % odd number
+            n=(n2+1)/2;
+            v=[b(1:n,:);b(n2:-1:n+1,:);NaN(1,32)];
+        else
+            n=n2/2;
+            v=[b(1:n,:);b(n2:-1:n+1,:)];
+        end
+        r=reshape(v,n,64);
+    end
+else
+    error('Not implemented for matrices')
 end
-
-
 % if isvector(a)
 %     n=numel(a)/32;
 %     b=reshape(a(:),2*n,16);
diff --git a/qem/qcheck.m b/qem/qcheck.m
index a4e8634d..06877db8 100644
--- a/qem/qcheck.m
+++ b/qem/qcheck.m
@@ -1,16 +1,26 @@
-function [dkk,stdf,meanf]=qcheck(k,k0,dang_ang,ax,largefig)
+function [dkk,qpdkk,dqdkk]=qcheck(qemres,k,k0,dang_ang,ax,largefig)
 %QCHECK		Analyzes quadrupole strength errors
 %
 %DKK=QCHECK(KL,KL0)
 %
+dqmask=qemres.atmodel.select('dq');
+quadmask=qemres.atmodel.select('qp');
+qpmask=false(length(qemres.at));
+qpmask(qemres.qpidx)=true;
+qpquadmask=qpmask & quadmask;
+qpdqmask=qpmask & dqmask;
 
 [nq,nm]=size(k);
-nfamq=nq;
 nfamd=length(dang_ang);
 labels=sr.qpname(1:17);
 
 k0l=k0(:,ones(1,nm));
 dkk=(k-k0l)./k0l;
+qpdkk=NaN(1,sum(quadmask));
+dqdkk=NaN(1,sum(dqmask));
+qpdkk(qpquadmask(quadmask))=dkk(qpquadmask(qpmask),end);
+dqdkk(qpdqmask(dqmask))=dkk(qpdqmask(qpmask),end);
+
 
 if  largefig
     
@@ -21,7 +31,7 @@ if  largefig
     %title(ax,'Gradient errors');
     ylabel(extfigax,'\DeltaK/K');
     grid(extfigax,'on');
-    set(extfigax,'Xlim',[0 nfamq+1],'FontSize',12);
+    set(extfigax,'Xlim',[0 nq+1],'FontSize',12);
     subplot(2,1,2);
     extfigax = gca;
     bar(extfigax,[dang_ang]);
@@ -32,27 +42,21 @@ if  largefig
     
 end
 
-Nc = length(dang_ang);
-extracor = mod(Nc,32);
-cor_per_cell_dip = squeeze(std(reshape(dang_ang((1+extracor/2):(Nc-extracor/2)),[],32),1,1));
-yyaxis(ax,'right');
-bar(ax,[cor_per_cell_dip],0.5,'EdgeColor','none');
-ylabel(ax,'\Delta\Theta/\Theta');
-grid(ax,'on');
-set(ax,'Xlim',[0 32+1],'FontSize',6);
-
-Nc = length(dkk);
-extracor = mod(Nc,32);
-cor_per_cell_quad = squeeze(std(reshape(dkk((1+extracor/2):(Nc-extracor/2)),[],32),1,1));
-yyaxis(ax,'left');
-bar(ax,[cor_per_cell_quad],0.9,'EdgeColor','none');
+% Nc = length(dang_ang);
+% extracor = mod(Nc,32);
+% cor_per_cell_dip = squeeze(std(reshape(dang_ang((1+extracor/2):(Nc-extracor/2)),[],32),1,1));
+% yyaxis(ax,'right');
+% bar(ax,[cor_per_cell_dip],0.5,'EdgeColor','none');
+% ylabel(ax,'\Delta\Theta/\Theta');
+% grid(ax,'on');
+% set(ax,'Xlim',[0 32+1],'FontSize',6);
+%cla(ax,'reset');
+stdquad=std(sr.fold(qpdkk),1,2,'omitnan');
+stddq=std(sr.fold(dqdkk),1,2,'omitnan');
+bar(ax,[stdquad(1:8);stddq],'EdgeColor','none');
 ylabel(ax,'\DeltaK/K');
 grid(ax,'on');
-set(ax,'Xlim',[0 32+1],'FontSize',6);
-ax.XTick=[1:32];
-ax.XTickLabel=arrayfun(@(a)num2str(a),[4:32,1:3],'un',0);
-ax.XTickLabelRotation=90;
-
+ax.XTickLabel={'QF1','QD2','QD3','QF4','QF4','QD5','QF6','QF8','DQ1','DQ2'};
 
 
 return
diff --git a/qem/qempaneldisp.m b/qem/qempaneldisp.m
index 8f1502e0..fc4d0a90 100644
--- a/qem/qempaneldisp.m
+++ b/qem/qempaneldisp.m
@@ -1,11 +1,18 @@
-function [qemb,qemres]=qempaneldisp(qemres,semres,qemb,qemb0,handles) %#ok<INUSL>
+function [qemb,qemres]=qempaneldisp(qemres,semres,qemb,qemb0,handles)
+
+cormask=false(length(qemres.at));
+cormask(qemres.qcoridx)=true;
+qpmask=false(length(qemres.at));
+qpmask(qemres.qpidx)=true;
+qpcormask=qpmask & cormask;
+
 tic
 if nargin < 5, handles=struct(); end
 
 if isfield(handles,'axes1')
-    qcheck(qemb.kn,qemb0.kn,qemb.dipdelta,handles.axes1,get(handles.largefigures,'Value'));
+    qcheck(qemres,qemb.kn,qemb0.kn,qemb.dipdelta,handles.axes1,get(handles.largefigures,'Value'));
 else
-    qcheck(qemb.kn,qemb0.kn,qemb.dipdelta,[],get(handles.largefigures,'Value'));
+    qcheck(qemres,qemb.kn,qemb0.kn,qemb.dipdelta,[],get(handles.largefigures,'Value'));
 end
 toc;
 qemb.at=qemat(qemres,qemb,true);    % builds the at structure
@@ -71,7 +78,7 @@ if all(isfield(qemres,{'resph','respv','frespx'}))
     if ~isempty(qemres.fithist) % store fit history
         qemres.fithist=[qemres.fithist;[std2(diffh(:)),std2(diffv(:)),hdispresidual]];
     else
-        qemres.fithist=[[std2(diffh(:)),std2(diffv(:)),hdispresidual]];
+        qemres.fithist=[std2(diffh(:)),std2(diffv(:)),hdispresidual];
     end
     
     fprintf('residual H = %g m/rad\n', std2(diffh(:)));
@@ -130,38 +137,45 @@ if isfield(handles,'axes4')                 % displays corr. strengths
         title(extfigax,'Corrector strengths');
         grid(extfigax,'on');
     end
-    Nc = length(qemb.cor);
-    extracor = mod(Nc,32); % remove extra correctors from injection cells
-    cor_per_cell_std = squeeze(std(reshape(qemb.cor((1+extracor/2):(Nc-extracor/2),:),[],32,size(qemb.cor,2)),1,1));
-    % cor_per_cell_mean = squeeze(mean(reshape(qemb.cor((1+extracor/2):(Nc-extracor/2),:),32,[],2),2));
-    wc= [0.9,0.5]; % superimpose ba
-    %set(gcf,'CurrentAxes',handles.axes4);
-    cla(handles.axes4);
-    hold off;
-    try
-        yyaxis rigth
-        bar(handles.axes4,cor_per_cell_std(:,1),wc(1),'EdgeColor','none');
-        
-        yyaxis left
-        bar(handles.axes4,cor_per_cell_std(:,2),wc(2),'EdgeColor','none');
-        yyaxis rigth
-    catch
-        bar(handles.axes4,cor_per_cell_std,0.9,'EdgeColor','none');    
+    scor=corstd(qemb.cor(:,1));
+    if size(qemb.cor,2) > 1
+        scor=[scor corstd(qemb.cor(:,end))];
     end
-    handles.axes4.XTick=[1:32];
-    handles.axes4.XTickLabel=arrayfun(@(a)num2str(a),[4:32,1:3],'un',0);
-    handles.axes4.XTickLabelRotation=90;
-    title(handles.axes4,'std Corrector strengths (per cell)');
+    bar(handles.axes4,scor(1:8,:));  % forget injection correctors
+        
+%     Nc = length(qemb.cor);
+%     extracor = mod(Nc,32); % remove extra correctors from injection cells
+%     cor_per_cell_std = squeeze(std(reshape(qemb.cor((1+extracor/2):(Nc-extracor/2),:),[],32,size(qemb.cor,2)),1,1));
+%     % cor_per_cell_mean = squeeze(mean(reshape(qemb.cor((1+extracor/2):(Nc-extracor/2),:),32,[],2),2));
+%     wc= [0.9,0.5]; % superimpose ba
+%     %set(gcf,'CurrentAxes',handles.axes4);
+%     cla(handles.axes4);
+%     hold off;
+%     try
+%         yyaxis rigth
+%         bar(handles.axes4,cor_per_cell_std(:,1),wc(1),'EdgeColor','none');
+%         
+%         yyaxis left
+%         bar(handles.axes4,cor_per_cell_std(:,2),wc(2),'EdgeColor','none');
+%         yyaxis rigth
+%     catch
+%         bar(handles.axes4,cor_per_cell_std,0.9,'EdgeColor','none');    
+%     end
+%     handles.axes4.XTick=[1:32];
+    handles.axes4.XTickLabel={'QF1','QD2','QD3','QF4','QF4','QD5','QF6','QF8'};
+%     handles.axes4.XTickLabelRotation=90;
+    ylabel(handles.axes4,'\Deltak/k');
+    title(handles.axes4,'std Corrector strengths (per family)');
     legend(handles.axes4,'initial','corrected');
     grid(handles.axes4,'on');
     
     
 end
 
-
-    function bave=baverage(b)
-        bb=mean(sr.fold(b),2);
-        bave=repmat([bb;bb(end:-1:1)],16,1);
+    function stdval=corstd(corval)
+        corstrength = corval./qemres.qcorl;
+        dkk=corstrength(qpcormask(cormask))./qemb0.kn(qpcormask(qpmask));
+        stdval=std(sr.fold(dkk),1,2,'omitnan');
     end
 
     function mdl=bmodul(beta,beta0)
diff --git a/qem/qempanelset.m b/qem/qempanelset.m
index a67001f9..190dd2de 100644
--- a/qem/qempanelset.m
+++ b/qem/qempanelset.m
@@ -26,7 +26,7 @@ catch %#ok<*CTCH>
     qemres=struct();
 end
 
-qemres.relative_or_absolute_fit = 'relative' ; 
+qemres.relative_or_absolute_fit = 'absolute' ; 
 
 qemres.bhscale=1;
 qemres.fithist = []; % history of fit
@@ -160,22 +160,29 @@ qemres.dipidx=find(dlmask | dqmask);
 % AnalyticResponseMatrix.m must be modified to include DQ! also qemderiv*.m
 % line 96 should be modified to get this working
 
-qpmask=atmod.select('qp');
+quadmask=atmod.select('qp');
 sfmask=atgetcells(atmod.ring,'FamName','S[FIJ]2[AE]');
 qf4mask=atgetcells(atmod.ring,'FamName','QF4[ABDEIJ]');
-fmask=qpmask | dqmask | sfmask;  % All QPs, DQs, SF2s
+% qpmask=quadmask;                   % All QPs
+% qpmask=quadmask | dqmask;          % All QPs, DQs
+qpmask=quadmask | dqmask | sfmask;	% All QPs, DQs, SF2s
 
-qemres.qpidx=find(fmask)';  % All QPs, DQs, SF2s
+qemres.qpidx=find(qpmask)';
 
 % % warning('only quadrupoles, no DQ!')
 qemres.quadidx = find(dqmask);	% Only quadrupoles
 
 qemres.dqidx=find(dqmask);
 
-qemres.qcoridx=find(qpmask)';    % Only quadrupoles
+cormask=quadmask;                   % Only quadrupoles
+qemres.qcoridx=find(cormask)';
 qemres.qcorl=atgetfieldvalues(qemres.at(qemres.qcoridx),'Length');
 
-qemres.qpfit = fmask(fmask) & ~qf4mask(fmask);
+qp_and_cor = qpmask & cormask;
+
+% qemres.qpfit = qpmask(qpmask);    % All true
+qemres.qpfit = ~sfmask(qpmask);     % No SF
+% qemres.qpfit = ~qf4mask(qpmask);  % No QF4
 qemres.ql=atgetfieldvalues(qemres.at(qemres.qpidx),'Length');
 
 qemres.sextidx=find(atmod.select('sx'));
@@ -251,7 +258,7 @@ if nargout >= 3
     qemb1.at=qemat(qemres,qemb1,false);
     
     if ComputeSimulatedRM
-        qemb(1)=qemoptics(qemres,qemb1);
+        qemb1=qemoptics(qemres,qemb1);
     else
         qemb1.lindata=NaN;       % BPM results
         qemb1.id07data=NaN;      % ID07 results
@@ -270,31 +277,34 @@ if nargout >= 3
         qemb1.frespx=NaN;
         qemb1.frespz=NaN;
         
-        qemb(1) = qemb1;
-        
     end
+    qemb(1) = qemb1;
+    
     
+    Nq=textscan(fopen(fullfile(qemres.datadir,'quadcor.dat'),'r'),'%s %f','HeaderLines',7);
+    indq=quadmask(quadmask | dqmask);
+    qemb(2).cor = Nq{2}(indq);	%select Quadrupoles
+%     Sq=textscan(fopen(fullfile(qemres.datadir,'skewcor.dat'),'r'),'%s %f','HeaderLines',7);
+%     qemb(2).skewcor = Sq{2};%zeros(length(qemres.skewidx),1);
+    qemb(2).skewcor=qemb1.skewcor;
     try
         [kn,qemb(2).dipdelta,qemb(2).dxs]=qemerrload(qemres,fullfile(qemres.opticsdir,'quadinierrors.mat'));
         qemb(2).kn=qemb(1).kn.*(1+srmodulation(kn));
     catch
         warning('qem:noinit','No initial errors in the optics directory');
-        qemb(2).kn=qemb(1).kn;
-        qemb(2).dipdelta=qemb(1).dipdelta;
-        qemb(2).dxs=qemb(1).dxs;
+        kn=qemb1.kn;
+        if strcmp(qemres.relative_or_absolute_fit,'absolute')
+            corval = qemb(2).cor./qemres.qcorl;
+            kn(qp_and_cor(qpmask)) = kn(qp_and_cor(qpmask)) - ...
+                corval(qp_and_cor(cormask));
+        end
+        qemb(2).dipdelta=qemb1.dipdelta;
+        qemb(2).dxs=qemb1.dxs;
     end
-    
-    Nq=textscan(fopen(fullfile(qemres.datadir,'quadcor.dat'),'r'),'%s %f','HeaderLines',7);
-    indq=qpmask(qpmask | dqmask);
-    corq = Nq{2}(indq); %select Quadrupoles
-    qemb(2).cor = corq;
-    Sq=textscan(fopen(fullfile(qemres.datadir,'skewcor.dat'),'r'),'%s %f','HeaderLines',7);
-    qemb(2).skewcor = Sq{2};%zeros(length(qemres.skewidx),1);
+    qemb(2).kn=kn;
     qemb(2).dzs=zeros(length(qemres.sextidx),1);
     qemb(2).diptilt=zeros(length(qemres.dipidx),1);
     qemb(2).ks=zeros(length(qemres.qpidx),1);
-    % mach=sr.opticsmatching(mach,'mux',76+qemres.tunes(1),'muy',27+qemres.tunes(2));      % retune the model
-    qemb(2).kn=atgetfieldvalues(mach(qemres.qpidx),'PolynomB',{2});
 end
 
 qemres=orderfields(qemres);
diff --git a/qem/qemresp.m b/qem/qemresp.m
index e5cbba37..7e3ee922 100644
--- a/qem/qemresp.m
+++ b/qem/qemresp.m
@@ -1,7 +1,6 @@
 function [rh,rv] = qemresp(mach,dct,bidx,sidx,elemfunc,varargin)
 %QEMRESP Builds the model response matrix
 %
-warning('off');
 nb=length(bidx);
 nq=length(sidx);
 vargs={[],[]};
diff --git a/sem/semtilt.m b/sem/semtilt.m
index 58519bdd..3df1a845 100644
--- a/sem/semtilt.m
+++ b/sem/semtilt.m
@@ -11,7 +11,7 @@ function [kl, tilt]=semtilt(klnorm,klskew)
 %
 oknorm=(klnorm~=0);
 kl=klskew;
-tilt=-pi/4*ones(size(klskew));
+tilt=-pi/4*(klskew~=0);
 ratio=klskew(oknorm)./klnorm(oknorm);
 kl(oknorm)=klnorm(oknorm).*sqrt(1+ratio.*ratio);
 tilt(oknorm)=-0.5*atan(ratio);
-- 
GitLab


From 8121fa1840b1ae60c852a066061cd6df63e569c0 Mon Sep 17 00:00:00 2001
From: Laurent Farvacque <laurent.farvacque@gmail.com>
Date: Mon, 6 Jan 2020 09:59:10 +0100
Subject: [PATCH 3/4] Bug fixes

---
 qem/qemguess.m     |  4 ++--
 qem/qempanel2.m    | 24 +++++++-----------------
 qem/qempanelgain.m |  2 +-
 qem/qempanelset.m  | 22 +++++++++++-----------
 qem/qemresp.m      | 10 +++++-----
 5 files changed, 26 insertions(+), 36 deletions(-)

diff --git a/qem/qemguess.m b/qem/qemguess.m
index a82ec4ab..4ba1fdd6 100644
--- a/qem/qemguess.m
+++ b/qem/qemguess.m
@@ -23,8 +23,8 @@ khg=khgain(qemres.hlist);
 kvr=kvrot(qemres.vlist);
 kvg=kvgain(qemres.vlist);
 
-goodh = all(~isnan(resph)); % not disabled horizontal steerers
-goodv = all(~isnan(respv)); % not disabled vertical steerers
+goodh = any(isfinite(resph)); % not disabled horizontal steerers
+goodv = any(isfinite(respv)); % not disabled vertical steerers
 
 bok=all(isfinite([...
     resph(:,goodh)...
diff --git a/qem/qempanel2.m b/qem/qempanel2.m
index 24bbf4af..b7385a9a 100644
--- a/qem/qempanel2.m
+++ b/qem/qempanel2.m
@@ -494,27 +494,17 @@ function pushbutton12_Callback(hObject, eventdata, handles)
 % handles    structure with handles and user data (see GUIDATA)
 global qemb qemres semres
 set(handles.statustext,'String','computing correction');
-fitproc={@qemcor_minus_fit,@qemcorrdt};
+fitproc={@qemcor_minus_fit,@qemcorrdt1};
 fitmethod=fitproc{get(handles.popupmenu1,'Value')};
 
-% get correctors starting point
-if strcmp(qemres.relative_or_absolute_fit,'absolute')
-        corval = qemb(2).cor(:,end);
-    elseif strcmp(qemres.relative_or_absolute_fit,'relative')
-        corval = qemb(2).cor(:,end)-qemb(2).cor(:,1);
-end
-
 if ~isempty(fitmethod)
-    
-    [newcor,qemres]=fitmethod(qemres,qemb(2).kn,corval);
-    set(handles.qemdata.cordef,'Enable','On');
-end
-
 % recompute total correction
-if strcmp(qemres.relative_or_absolute_fit,'absolute')
-    qemb(2).cor(:,2) = newcor;
-elseif strcmp(qemres.relative_or_absolute_fit,'relative')
-    qemb(2).cor(:,2) = qemb(2).cor(:,1) + newcor;
+    if strcmp(qemres.relative_or_absolute_fit,'absolute')
+        [qemb(2).cor(:,2),qemres] = fitmethod(qemres,qemb(2).kn);
+    elseif strcmp(qemres.relative_or_absolute_fit,'relative')
+        [qemb(2).cor(:,2),qemres] = fitmethod(qemres,qemb(2).kn, qemb(2).cor(:,end));
+    end
+    set(handles.qemdata.cordef,'Enable','On');
 end
 
 [qemb(2),qemres]=qempaneldisp(qemres,semres,qemb(2),qemb(1),handles);
diff --git a/qem/qempanelgain.m b/qem/qempanelgain.m
index 77df8f8c..371f8cb3 100644
--- a/qem/qempanelgain.m
+++ b/qem/qempanelgain.m
@@ -1,7 +1,7 @@
 function [kgain,lgain]=qempanelgain(measresp,modelresp,kgini,bgini)
 
 ok=isfinite(measresp);
-goodsteerer = all(~isnan(measresp) & ~isinf(measresp));
+goodsteerer = any(isfinite(measresp));
 bpmlist=sum(ok,2)>size(measresp(:,goodsteerer,1),2)-1; % sum only columns that are not NaN (disabled/frozen/fault correctors) 
 %nok=224-sum(bgini(~bpmlist));
 sumgain=sum(bgini(bpmlist));
diff --git a/qem/qempanelset.m b/qem/qempanelset.m
index 190dd2de..a895ab5f 100644
--- a/qem/qempanelset.m
+++ b/qem/qempanelset.m
@@ -22,7 +22,7 @@ end
 initname=fullfile(MACHFS,'MDT','respmat','qeminit.mat');
 try
     qemres=load(initname);
-catch %#ok<*CTCH>
+catch
     qemres=struct();
 end
 
@@ -37,7 +37,6 @@ qemres.opticsdir=fullfile(dirname,'optics','');
 try
     simu = tango.Device('sys/ringsimulator/ebs');
     qemres.simufile=simu.get_property('RingFile');
-    
 catch
     disp('tango ''sys/ringsimulator/ebs'' not available')
     qemres.simufile ='None';
@@ -45,7 +44,7 @@ end
 
 qemres.opticsname = theodir;
 a=load(fullfile(qemres.opticsdir,'betamodel.mat'),'betamodel');
-ring= a.betamodel;
+ring=atradoff(a.betamodel,'IdentityPass','auto','auto');
 
 warning('check SB are here')
 
@@ -165,12 +164,12 @@ sfmask=atgetcells(atmod.ring,'FamName','S[FIJ]2[AE]');
 qf4mask=atgetcells(atmod.ring,'FamName','QF4[ABDEIJ]');
 % qpmask=quadmask;                   % All QPs
 % qpmask=quadmask | dqmask;          % All QPs, DQs
-qpmask=quadmask | dqmask | sfmask;	% All QPs, DQs, SF2s
+errmask=quadmask | dqmask | sfmask;	% All QPs, DQs, SF2s
 
-qemres.qpidx=find(qpmask)';
+qemres.qpidx=find(errmask)';
 
 % % warning('only quadrupoles, no DQ!')
-qemres.quadidx = find(dqmask);	% Only quadrupoles
+qemres.quadidx = find(quadmask);	% Only quadrupoles
 
 qemres.dqidx=find(dqmask);
 
@@ -178,11 +177,12 @@ cormask=quadmask;                   % Only quadrupoles
 qemres.qcoridx=find(cormask)';
 qemres.qcorl=atgetfieldvalues(qemres.at(qemres.qcoridx),'Length');
 
-qp_and_cor = qpmask & cormask;
+qp_and_cor = errmask & cormask;
 
-% qemres.qpfit = qpmask(qpmask);    % All true
-qemres.qpfit = ~sfmask(qpmask);     % No SF
-% qemres.qpfit = ~qf4mask(qpmask);  % No QF4
+% qemres.qpfit = qpmask(errmask);         % All true
+% qemres.qpfit = ~sfmask(errmask);        % No SF
+% qemres.qpfit = ~(sfmask(errmask) | dqmask(errmask));	  % No DQ, no SF
+qemres.qpfit = ~qf4mask(errmask);       % No QF4
 qemres.ql=atgetfieldvalues(qemres.at(qemres.qpidx),'Length');
 
 qemres.sextidx=find(atmod.select('sx'));
@@ -295,7 +295,7 @@ if nargout >= 3
         kn=qemb1.kn;
         if strcmp(qemres.relative_or_absolute_fit,'absolute')
             corval = qemb(2).cor./qemres.qcorl;
-            kn(qp_and_cor(qpmask)) = kn(qp_and_cor(qpmask)) - ...
+            kn(qp_and_cor(errmask)) = kn(qp_and_cor(errmask)) - ...
                 corval(qp_and_cor(cormask));
         end
         qemb(2).dipdelta=qemb1.dipdelta;
diff --git a/qem/qemresp.m b/qem/qemresp.m
index 7e3ee922..e60c587b 100644
--- a/qem/qemresp.m
+++ b/qem/qemresp.m
@@ -7,7 +7,7 @@ vargs={[],[]};
 vargs(1:length(varargin))=varargin;
 orbit0=vargs{1};
 kick=0.000025*ones(1,nq);       % default kick: 25 urad
-if isempty(orbit0), orbit0=findsyncorbit(mach,dct,bidx); end
+%if isempty(orbit0), orbit0=findsyncorbit(mach,dct,bidx); end
 if ~isempty(vargs{2}), kick(:)=vargs{2}; end
 rh=zeros(nb,nq);
 rv=zeros(nb,nq);
@@ -17,11 +17,11 @@ for iq=1:nq
     mach{idq}=elemfunc(esave,iq,kick(iq));
     orbitp=findsyncorbit(mach,dct,bidx);
     %orbitp=findorbit6Err(mach,bidx);
-    %mach{idq}=elemfunc(esave,iq,-kick(iq)); %double sided computation
-    %orbitm=findsyncorbit(mach,dct,bidx);
+    mach{idq}=elemfunc(esave,iq,-kick(iq)); %double sided computation
+    orbitm=findsyncorbit(mach,dct,bidx);
     %orbitm=findorbit6Err(mach,bidx);
-    %orbit=(orbitp-orbitm)/kick(iq)/2;
-    orbit=(orbitp-orbit0)/kick(iq);
+    orbit=(orbitp-orbitm)/kick(iq)/2;
+    %orbit=(orbitp-orbit0)/kick(iq);
     mach{idq}=esave;
     rh(:,iq)=orbit(1,:)';
     rv(:,iq)=orbit(3,:)';
-- 
GitLab


From 0abf91443466b67d477deb640d0c58489073b234 Mon Sep 17 00:00:00 2001
From: Laurent Farvacque <laurent.farvacque@esrf.fr>
Date: Tue, 7 Jan 2020 15:33:10 +0100
Subject: [PATCH 4/4] merged from master

---
 qem/qemderivAnalytic.m |  2 +-
 qem/qempanel2.m        |  2 +-
 qem/qempanelset.m      | 13 +++++++++++++
 3 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/qem/qemderivAnalytic.m b/qem/qemderivAnalytic.m
index d6a5c7ac..58e22c46 100644
--- a/qem/qemderivAnalytic.m
+++ b/qem/qemderivAnalytic.m
@@ -88,7 +88,7 @@ qind=find(atgetcells(mach,'Class','Quadrupole') )'; % | atgetcells(mach,'FamName
 % if isfield(qemres,'quadhdispresponse')
 %     % do not recomput if already available
 %     dDx_dq=qemres.quadhdispresponse;
-%     
+%
 % else
     disp('FIRST and only h.disp derivative with respect to quadrupoles');
     disp('computation starts on lattice without errors/correctors');
diff --git a/qem/qempanel2.m b/qem/qempanel2.m
index b7385a9a..d60dd2a5 100644
--- a/qem/qempanel2.m
+++ b/qem/qempanel2.m
@@ -304,7 +304,7 @@ if get(handles.DCRMcheckbox,'Value') % DC/AC DC measurement
         % if full measurement RM CSV files will be stored. also disabled magnets are in the list, an empty response is stored.
         rc.MeasureOrbitResponseMatrix('full',true); 
     else
-        rc.MeasureOrbitResponseMatrix('Hsteerer',steerlist,'Vsteerer',qemres.vlist); 
+        rc.MeasureOrbitResponseMatrix('Hsteerer',steerlist,'Vsteerer',qemres.vlist);
     end
     
     % store measured RM in theory folder.
diff --git a/qem/qempanelset.m b/qem/qempanelset.m
index a895ab5f..2e56ffe5 100644
--- a/qem/qempanelset.m
+++ b/qem/qempanelset.m
@@ -45,6 +45,19 @@ end
 qemres.opticsname = theodir;
 a=load(fullfile(qemres.opticsdir,'betamodel.mat'),'betamodel');
 ring=atradoff(a.betamodel,'IdentityPass','auto','auto');
+% warning('setting inj bump')
+% warning('setting inj bump')
+% warning('setting inj bump')
+% warning('setting inj bump')
+% a = load('/machfs/MDT/2019/2019_12_16/str.mat');
+% indVst = find(atgetcells(ring,'FamName','S[HDFIJ]\w*'))';
+%
+% ring = atsetfieldvalues(ring,indVst([1:4,end-3:end]),'KickAngle',{1,1},a.str([1:4,end-3:end])');
+% warning('setting inj bump')
+% warning('setting inj bump')
+% warning('setting inj bump')
+% warning('setting inj bump')
+% warning('setting inj bump')
 
 warning('check SB are here')
 
-- 
GitLab