From 5d402eb15e0e027b59d71ea29a78f61d513deac9 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Mon, 26 Jun 2023 10:13:10 +0200 Subject: [PATCH 01/26] When solving system of linear equations to get beam size or position, no longer use linsolve but lsqminnorm, such that (in case of non-unique solution) the solution with min module is chosen --- optics/SolveOrbSystem.m | 19 +++++++++---------- optics/SolveSigSystem.m | 20 +++++++++----------- 2 files changed, 18 insertions(+), 21 deletions(-) diff --git a/optics/SolveOrbSystem.m b/optics/SolveOrbSystem.m index d6c0201..a098267 100644 --- a/optics/SolveOrbSystem.m +++ b/optics/SolveOrbSystem.m @@ -57,15 +57,14 @@ end function X=SolveOrbSystemActual(A,BARsM,lDebug) - opts.RECT=true; - [X,r]=linsolve(A,BARsM,opts); - if ( size(A,2)==3 && r==2 && sum(A(:,3))>0 ) - if ( lDebug ) - fprintf("==> new ORB fit %gD, %g points\n",size(A,2),size(A,1)); - end - X=linsolve(A,BARsM,opts); % let MatLab raise the warning - if ( lDebug ) - fprintf("==> done.\n"); - end + if ( lDebug ) + fprintf("==> new ORB fit %gD, %g points\n",size(A,2),size(A,1)); + end + % solves the linear equation AX = B and minimizes the value of norm(A*X-B). + % If several solutions exist to this problem, then lsqminnorm returns the + % solution that minimizes norm(X) + X=lsqminnorm(A,BARsM); + if ( lDebug ) + fprintf("==> done.\n"); end end diff --git a/optics/SolveSigSystem.m b/optics/SolveSigSystem.m index 32c00c6..1f2a55c 100644 --- a/optics/SolveSigSystem.m +++ b/optics/SolveSigSystem.m @@ -84,18 +84,16 @@ end function X=SolveSigSystemActual(A,sigs2,initConds,lDebug) - opts.RECT=true; if ( ismissing(initConds) ) - % X=linsolve(A,sigs2,opts); - [X,r]=linsolve(A,sigs2,opts); - if ( ( size(A,2)==5 && r==3 && sum(A(:,4:5),"all")>0 ) || ( size(A,2)==3 && r==2 ) ) - if ( lDebug ) - fprintf("==> new SIG fit %gD, %g points\n",size(A,2),size(A,1)); - end - X=linsolve(A,sigs2,opts); % let MatLab raise the warning - if ( lDebug ) - fprintf("==> done.\n"); - end + if ( lDebug ) + fprintf("==> new SIG fit %gD, %g points\n",size(A,2),size(A,1)); + end + % solves the linear equation AX = B and minimizes the value of norm(A*X-B). + % If several solutions exist to this problem, then lsqminnorm returns the + % solution that minimizes norm(X) + X=lsqminnorm(A,sigs2); % let MatLab raise the warning + if ( lDebug ) + fprintf("==> done.\n"); end else % lb=[0 -inf 0]; ub=[inf inf inf]; XC = lsqlin(A,sigs2,[],[],[],[],lb,ub); From 405160d1244631d9f36abd1aab472d7d1fc86629 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Mon, 26 Jun 2023 10:16:10 +0200 Subject: [PATCH 02/26] adding possibility to parse .tfs with scan currents --- MADX-optics/GetColumnsAndMappingTFS.m | 214 +++++++++++++++----------- MADX-optics/ParseTfsTable.m | 44 ++++-- 2 files changed, 158 insertions(+), 100 deletions(-) diff --git a/MADX-optics/GetColumnsAndMappingTFS.m b/MADX-optics/GetColumnsAndMappingTFS.m index 003da87..7a913b2 100644 --- a/MADX-optics/GetColumnsAndMappingTFS.m +++ b/MADX-optics/GetColumnsAndMappingTFS.m @@ -1,19 +1,23 @@ function [ colNames, colUnits, colFacts, mapping, readFormat ] = ... - GetColumnsAndMappingTFS(whichTFS,genBy) + GetColumnsAndMappingTFS(whichTFS,genBy,fileName) % GetColumnsAndMappingTFS return column names, units and column -% mapping of twiss tables +% mapping of (predefined) twiss tables % % [ colNames, colUnits, colFacts, mapping, readFormat ] = % GetColumnsAndMappingTFS(whichTFS) % % input arguments: -% whichData: format of the table: +% whichTFS: format of the table: % 'OPTICS': optics table by MAD-X; % 'GEOMETRY': lattice geometry by MAD-X; % 'RMATRIX': (selected elements of) response matrix by MAD-X; +% 'CURR': a TFS table listing scan currents; % genBy: format of the table: % 'TWISS': a TFS table as generated by a MADX TWISS command; % 'SCAN': a TFS table as generated by a current scan; +% fileName (optional): name of file; +% the file is parsed for auto-detection of format. +% FOR THE TIME BEING, ONLY "CURR" TYPE IS SUPPORTED! % % output arguments: % colNames: array with name of columns/variables; @@ -23,7 +27,7 @@ % readFormat: format string necessary to correctly parse the file; % % This function states that the optics, geometry and Rmatrix TFS tables -% have specific formats: +% have specific pre-defined formats: % - generated by a MADX TWISS command (genBy=="TWISS") % . optics: % NAME, KEYWORD, L, S, BETX, ALFX, BETY, ALFY, X, PX, Y, PY, @@ -34,7 +38,6 @@ % . response matrix: % NAME, KEYWORD, L, S, RE11, RE12, RE21, RE22, RE16, RE26, RE33, RE34, % RE43, RE44, RE36, RE46, RE51, RE52, RE55, RE56, RE66; -% % - generated by a scan (genBy=="SCAN") % . optics: % Brho[Tm], BP[mm], myID[], BETX, ALFX, BETY, ALFY, X, PX, Y, PY, @@ -47,90 +50,129 @@ % RE43, RE44, RE36, RE46, RE51, RE52, RE55, RE56, RE66; % % See also ParseTfsTable. - if (~exist("genBy","var") | ismissing(genBy)), genBy="TWISS"; end % default genBy + % default genBy + if (~exist("genBy","var") | ismissing(genBy)) + if (strcmpi(whichTFS,"CURR")) + genBy="SCAN"; + else + genBy="TWISS"; + end + end % from label to data column - switch upper(genBy) - case "TWISS" - switch upper(whichTFS) - case {'OPT','OPTICS'} - colNames=[ "NAME" "KEYWORD" "L" "S" "BETX" "ALFX" "BETY" "ALFY" ... - "X" "PX" "Y" "PY" "DX" "DPX" "DY" "DPY" ... - "MUX" "MUY" ]; - colUnits=[ "" "" "m" "m" "m" "" "m" "" ... - "m" "" "m" "" "m" "" "m" "" ... - "2\pi" "2\pi"]; - readFormat = '%s %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f'; - case {'GEO','GEOMETRY'} - colNames=[ "NAME" "KEYWORD" "L" "S" ... - "KICK" "HKICK" "VKICK" "ANGLE" "K0L" "K1L" "K2L" ... - "APERTYPE" "APER_1" "APER_2" "APER_3" "APER_4" "APOFF_1" "APOFF_2" ]; - colUnits=[ "" "" "m" "m" ... - "rad" "rad" "rad" "rad" "rad" "rad m^{-1}" "rad m^{-2}" ... - "" "m" "m" "m" "m" "m" "m" ]; - readFormat = '%s %s %f %f %f %f %f %f %f %f %f %s %f %f %f %f %f %f %f'; - case {'RM','RMATRIX'} - colNames=[ "NAME" "KEYWORD" "L" "S" ... - "RE11" "RE12" "RE21" "RE22" "RE16" "RE26" ... - "RE33" "RE34" "RE43" "RE44" "RE36" "RE46" ... - "RE51" "RE52" "RE55" "RE56" "RE66" ]; - colUnits=[ "" "" "m" "m" ... - "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... - "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... - "s m^{-1}" "s rad^{-1}" "" "s" "m" ]; - readFormat = '%s %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f'; - otherwise - error('which column mapping for TFS table?'); - end - case "SCAN" - mySep=','; - switch upper(whichTFS) - case {'OPT','OPTICS'} - colNames=[ "BRHO" "BP" "ID" ... - "BETX" "ALFX" "BETY" "ALFY" ... - "X" "PX" "Y" "PY" ... - "DX" "DPX" "DY" "DPY" ... - "MUX" "MUY" ]; - colUnits=[ "Tm" "mm" "" ... - "m" "" "m" "" ... - "m" "" "m" "" ... - "m" "" "m" "" ... - "2\pi" "2\pi"]; - readFormat = strcat('%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f' ); - case {'GEO','GEOMETRY'} - colNames=[ "BRHO" "BP" "ID" ... - "KICK" "HKICK" "VKICK" "ANGLE" "K0L" "K1L" "K2L" ... - "APERTYPE" "APER_1" "APER_2" "APER_3" "APER_4" "APOFF_1" "APOFF_2" ]; - colUnits=[ "Tm" "mm" "" ... - "rad" "rad" "rad" "rad" "rad" "rad m^{-1}" "rad m^{-2}" ... - "" "m" "m" "m" "m" "m" "m" ]; - readFormat = strcat('%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f' ); - case {'RM','RMATRIX'} - colNames=[ "BRHO" "BP" "ID" ... - "RE11" "RE12" "RE21" "RE22" "RE16" "RE26" ... - "RE33" "RE34" "RE43" "RE44" "RE36" "RE46" ... - "RE51" "RE52" "RE55" "RE56" "RE66" ]; - colUnits=[ "Tm" "mm" "" ... - "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... - "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... - "s m^{-1}" "s rad^{-1}" "" "s" "m" ]; - readFormat = strcat('%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... - '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f' ); - otherwise - error('which column mapping for TFS table?'); + if (exist("fileName","var")) + fprintf("auto-detecting header of file %s ...\n",fileName); + if (~strcmpi(whichTFS,"CURR")) + error("auto-detection of file header only available for the CURR type!"); + end + fileID = fopen(fileName,'r'); + colNames=missing(); + while ~feof(fileID) + tline=strip(fgetl(fileID)); % drop new-line char + if (startsWith(tline,"*")) + colNames=split(tline); + colNames=colNames(strlength(colNames)>0); % forget about empty strings due to consecutive blanks + colNames=colNames(2:end); % forget about heading "*" + colNames=strrep(colNames,"_A",""); % drop closing _A + break end - otherwise - error("How was the .tfs table generated (TWISS/SCAN)? %s NOT recognised!",genBy); + end + fclose(fileID); + if (ismissing(colNames)) + error("No line headed by '*' char was found. Are you sure about the file format?"); + end + colUnits=strings(length(colNames),1); + colUnits(:)="A"; + readFormat=strings(1,length(colNames)); + readFormat(:)="%f"; + readFormat=join(readFormat); + else + % pre-defined formats + if (strcmpi(whichTFS,"CURR")) + error("NO pre-defined format for CURR .tfs files!"); + end + switch upper(genBy) + case "TWISS" + switch upper(whichTFS) + case {'OPT','OPTICS'} + colNames=[ "NAME" "KEYWORD" "L" "S" "BETX" "ALFX" "BETY" "ALFY" ... + "X" "PX" "Y" "PY" "DX" "DPX" "DY" "DPY" ... + "MUX" "MUY" ]; + colUnits=[ "" "" "m" "m" "m" "" "m" "" ... + "m" "" "m" "" "m" "" "m" "" ... + "2\pi" "2\pi"]; + readFormat = '%s %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f'; + case {'GEO','GEOMETRY'} + colNames=[ "NAME" "KEYWORD" "L" "S" ... + "KICK" "HKICK" "VKICK" "ANGLE" "K0L" "K1L" "K2L" ... + "APERTYPE" "APER_1" "APER_2" "APER_3" "APER_4" "APOFF_1" "APOFF_2" ]; + colUnits=[ "" "" "m" "m" ... + "rad" "rad" "rad" "rad" "rad" "rad m^{-1}" "rad m^{-2}" ... + "" "m" "m" "m" "m" "m" "m" ]; + readFormat = '%s %s %f %f %f %f %f %f %f %f %f %s %f %f %f %f %f %f %f'; + case {'RM','RMATRIX'} + colNames=[ "NAME" "KEYWORD" "L" "S" ... + "RE11" "RE12" "RE21" "RE22" "RE16" "RE26" ... + "RE33" "RE34" "RE43" "RE44" "RE36" "RE46" ... + "RE51" "RE52" "RE55" "RE56" "RE66" ]; + colUnits=[ "" "" "m" "m" ... + "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... + "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... + "s m^{-1}" "s rad^{-1}" "" "s" "m" ]; + readFormat = '%s %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f'; + otherwise + error('which column mapping for TFS table?'); + end + case "SCAN" + mySep=','; + switch upper(whichTFS) + case {'OPT','OPTICS'} + colNames=[ "BRHO" "BP" "ID" ... + "BETX" "ALFX" "BETY" "ALFY" ... + "X" "PX" "Y" "PY" ... + "DX" "DPX" "DY" "DPY" ... + "MUX" "MUY" ]; + colUnits=[ "Tm" "mm" "" ... + "m" "" "m" "" ... + "m" "" "m" "" ... + "m" "" "m" "" ... + "2\pi" "2\pi"]; + readFormat = strcat('%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f' ); + case {'GEO','GEOMETRY'} + colNames=[ "BRHO" "BP" "ID" ... + "KICK" "HKICK" "VKICK" "ANGLE" "K0L" "K1L" "K2L" ... + "APERTYPE" "APER_1" "APER_2" "APER_3" "APER_4" "APOFF_1" "APOFF_2" ]; + colUnits=[ "Tm" "mm" "" ... + "rad" "rad" "rad" "rad" "rad" "rad m^{-1}" "rad m^{-2}" ... + "" "m" "m" "m" "m" "m" "m" ]; + readFormat = strcat('%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f' ); + case {'RM','RMATRIX'} + colNames=[ "BRHO" "BP" "ID" ... + "RE11" "RE12" "RE21" "RE22" "RE16" "RE26" ... + "RE33" "RE34" "RE43" "RE44" "RE36" "RE46" ... + "RE51" "RE52" "RE55" "RE56" "RE66" ]; + colUnits=[ "Tm" "mm" "" ... + "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... + "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... + "s m^{-1}" "s rad^{-1}" "" "s" "m" ]; + readFormat = strcat('%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... + '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f' ); + otherwise + error('which column mapping for TFS table?'); + end + otherwise + error("How was the .tfs table generated (TWISS/SCAN)? %s NOT recognised!",genBy); + end end - colFacts = ones(1,length(colNames)); - mapping = 1:length(colNames); + colFacts = ones(1,length(colNames)); + mapping = 1:length(colNames); end diff --git a/MADX-optics/ParseTfsTable.m b/MADX-optics/ParseTfsTable.m index 67d6578..9015985 100644 --- a/MADX-optics/ParseTfsTable.m +++ b/MADX-optics/ParseTfsTable.m @@ -23,34 +23,50 @@ % genBy: format of the table: % 'TWISS': a TFS table as generated by a MADX TWISS command; % 'SCAN': a TFS table as generated by a current scan; +% 'CURR': a TFS table listing scan currents; % Nota Bene: the column names and formats are considered as part of the % header; % % See also GetColumnsAndMappingTFS, ParseOpticsFileHeader, GetAperture, PlotLattice. - if (~exist("genBy","var") | ismissing(genBy)), genBy="TWISS"; end % default genBy - switch upper(genBy) - case "TWISS" - nHeader=48; - case "SCAN" - nHeader=1; - otherwise - error("How was the .tfs table generated (TWISS/SCAN)? %s NOT recognised!",genBy); + if (~exist("genBy","var") | ismissing(genBy)) + if (strcmpi(whichTFS,"CURR")) + genBy="SCAN"; + else + genBy="TWISS"; + end + end + if (strcmpi(whichTFS,"CURR")) + nHeader=7; + else + switch upper(genBy) + case "TWISS" + nHeader=48; + case "SCAN" + nHeader=1; + otherwise + error("How was the .tfs table generated (TWISS/SCAN)? %s NOT recognised!",genBy); + end end - - % get format of the file - [ colNames, colUnits, colFacts, mapping, readFormat ] = ... - GetColumnsAndMappingTFS(whichTFS,genBy); if ( length(fileNames)==1 ) fprintf('parsing %s file %s generated by %s ...\n',whichTFS,fileNames,genBy); + % - get format of the file + if (strcmpi(whichTFS,"CURR")) + [ colNames, colUnits, colFacts, mapping, readFormat ] = ... + GetColumnsAndMappingTFS(whichTFS,genBy,fileNames); + else + [ colNames, colUnits, colFacts, mapping, readFormat ] = ... + GetColumnsAndMappingTFS(whichTFS,genBy); + end + % - actually parse file fileID = fopen(fileNames,'r'); tfsTable = textscan(fileID,readFormat,'HeaderLines',nHeader); fclose(fileID); else tfsTable=[]; - for fileName=fileNames - tmpTfsTable = ParseTfsTable(fileName,whichTFS,genBy); + for iFileName=1:length(fileNames) + tmpTfsTable = ParseTfsTable(fileNames(iFileName),whichTFS,genBy); tfsTable=[ tfsTable ; tmpTfsTable ]; end fprintf('...acquired %i %s files.\n',length(fileNames),whichTFS); From 318c3466c51bb9bfbf5a6c2aa723e8d81b9d76b6 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Mon, 26 Jun 2023 10:17:31 +0200 Subject: [PATCH 03/26] adding minor utilities --- MADX-optics/ShowOpticsTFSatEle.m | 110 +++++++++++++++++++++++++++++++ general/BinCentres2Edges.m | 10 +++ operations/QAxls2tfs.m | 23 +++++++ 3 files changed, 143 insertions(+) create mode 100644 MADX-optics/ShowOpticsTFSatEle.m create mode 100644 general/BinCentres2Edges.m create mode 100644 operations/QAxls2tfs.m diff --git a/MADX-optics/ShowOpticsTFSatEle.m b/MADX-optics/ShowOpticsTFSatEle.m new file mode 100644 index 0000000..0287a7b --- /dev/null +++ b/MADX-optics/ShowOpticsTFSatEle.m @@ -0,0 +1,110 @@ +function ShowOpticsTFSatEle(myOptics,optLabels) + %% identify columns + [ colNames, colUnits, colFacts, mapping, readFormat ] = ... + GetColumnsAndMappingTFS("optics","SCAN"); + planes=["x" "y"]; + + nPoints=length(myOptics{1}); + nOptics=size(myOptics,1); + if (nOptics==1) + nRows=2; nCols=2; + else + nRows=4; nCols=2; + if (~exist("optLabels","var")) + optLabels=compose("optics #%d",1:nOptics); + end + end + + %% actually plot + iPlot=0; + figure(); + cm=colormap(jet(nOptics)); + iCol=[ mapping(strcmp(colNames,'BETX')) mapping(strcmp(colNames,'BETY')) ]; + if (nOptics==1) + iPlot=iPlot+1; + subplot(nRows,nCols,iPlot); + plot(1:nPoints,myOptics{nOptics,iCol(1)},"ro-",1:nPoints,myOptics{nOptics,iCol(2)},"b*-"); + ylabel("\beta [m]"); grid on; xlabel("ID []"); + legend(compose("\\beta_{%s}",planes),"Location","best"); + else + for iPlane=1:2 + % planes + iPlot=iPlot+1; + subplot(nRows,nCols,iPlot); + for iOpt=1:nOptics + % optics + if (iOpt>1), hold on; end + plot(1:nPoints,myOptics{iOpt,iCol(iPlane)},".-","Color",cm(iOpt,:)); + end + ylabel(sprintf("\\beta_{%s} [m]",planes(iPlane))); grid on; xlabel("ID []"); + legend(optLabels,"Location","best"); + end + end + + iCol=[ mapping(strcmp(colNames,'X')) mapping(strcmp(colNames,'Y')) ]; + if (nOptics==1) + iPlot=iPlot+1; + subplot(nRows,nCols,iPlot); + plot(1:nPoints,myOptics{nOptics,iCol(1)},"ro-",1:nPoints,myOptics{nOptics,iCol(2)},"b*-"); + ylabel("pos [m]"); grid on; xlabel("ID []"); + legend(compose("%s",planes),"Location","best"); + else + for iPlane=1:2 + % planes + iPlot=iPlot+1; + subplot(nRows,nCols,iPlot); + for iOpt=1:nOptics + % optics + if (iOpt>1), hold on; end + plot(1:nPoints,myOptics{iOpt,iCol(iPlane)},".-","Color",cm(iOpt,:)); + end + ylabel(sprintf("%s [m]",planes(iPlane))); grid on; xlabel("ID []"); + legend(optLabels,"Location","best"); + end + end + + iCol=[ mapping(strcmp(colNames,'DX')) mapping(strcmp(colNames,'DY')) ]; + if (nOptics==1) + iPlot=iPlot+1; + subplot(nRows,nCols,iPlot); + plot(1:nPoints,myOptics{nOptics,iCol(1)},"ro-",1:nPoints,myOptics{nOptics,iCol(2)},"b*-"); + ylabel("D [m]"); grid on; xlabel("ID []"); + legend(compose("D_{%s}",planes),"Location","best"); + else + for iPlane=1:2 + % planes + iPlot=iPlot+1; + subplot(nRows,nCols,iPlot); + for iOpt=1:nOptics + % optics + if (iOpt>1), hold on; end + plot(1:nPoints,myOptics{iOpt,iCol(iPlane)},".-","Color",cm(iOpt,:)); + end + ylabel(sprintf("D_{%s} [m]",planes(iPlane))); grid on; xlabel("ID []"); + legend(optLabels,"Location","best"); + end + end + + iCol=[ mapping(strcmp(colNames,'MUX')) mapping(strcmp(colNames,'MUY')) ]; + if (nOptics==1) + iPlot=iPlot+1; + subplot(nRows,nCols,iPlot); + plot(1:nPoints,myOptics{nOptics,iCol(1)}*360,"ro-",1:nPoints,myOptics{nOptics,iCol(2)}*360,"b*-"); + ylabel("\mu [degs]"); grid on; xlabel("ID []"); + legend(compose("\\mu_{%s}",planes),"Location","best"); + else + for iPlane=1:2 + % planes + iPlot=iPlot+1; + subplot(nRows,nCols,iPlot); + for iOpt=1:nOptics + % optics + if (iOpt>1), hold on; end + plot(1:nPoints,myOptics{iOpt,iCol(iPlane)}*360,".-","Color",cm(iOpt,:)); + end + ylabel(sprintf("\\mu_{%s} [degs]",planes(iPlane))); grid on; xlabel("ID []"); + legend(optLabels,"Location","best"); + end + end + +end diff --git a/general/BinCentres2Edges.m b/general/BinCentres2Edges.m new file mode 100644 index 0000000..d1fbb65 --- /dev/null +++ b/general/BinCentres2Edges.m @@ -0,0 +1,10 @@ +function [Vout]=BinCentres2Edges(Vin) + delta=diff(Vin); delta=delta(1); + if (size(Vin,1)==1) + % column array + Vout=[Vin Vin(end)+delta]-delta/2; + else + % row array + Vout=[Vin;Vin(end)+delta]-delta/2; + end +end diff --git a/operations/QAxls2tfs.m b/operations/QAxls2tfs.m new file mode 100644 index 0000000..adfe236 --- /dev/null +++ b/operations/QAxls2tfs.m @@ -0,0 +1,23 @@ +function QAxls2tfs(iFileName,oFileName) + if (~exist("oFileName","var")) + [tmpDir,tmpNam,tmpExt]=fileparts(iFileName); + oFileName=strcat(tmpDir,"/",tmpNam,".tfs"); + end + fprintf("Converting QA file %s into TFS-like file %s ...\n",iFileName,oFileName); + + %% get data and store it + myData=readcell(iFileName); + PSNames=string(myData(2:end,1)); + currValues=cell2mat(myData(2:end,3:end)); + + %% convert to .tfs format + headers=compose("%s_A",PSNames); + headerTypes=strings(1,length(headers)); + headerTypes(:)="%le"; + myTitle=sprintf("from %s",iFileName); + myTable=currValues'; + ExportMADXtable(oFileName,myTitle,myTable,headers,headerTypes); + + %% + fprintf("...done;\n"); +end From a8a1874e584cd6da33310f8824a46bb6226c3a2b Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 5 Jul 2023 14:14:35 +0200 Subject: [PATCH 04/26] Possibility to change on the fly what should be displaied as independent variable in scans/profiles --- DisplayBeamProfiles.m | 33 ++++++++++++++----- ...tingConditions.m => StartConditionsRead.m} | 0 2 files changed, 25 insertions(+), 8 deletions(-) rename MADX-tracking/{ReadStartingConditions.m => StartConditionsRead.m} (100%) diff --git a/DisplayBeamProfiles.m b/DisplayBeamProfiles.m index bd7f908..5e3760a 100644 --- a/DisplayBeamProfiles.m +++ b/DisplayBeamProfiles.m @@ -53,6 +53,7 @@ MonPaths=[... strcat(kPath,"\Area dati MD\00Summary\Carbonio\2023\Maggio\2023.05.09-10\Steering ridotti\GIM\PRC-544-230511-0028_H2-009B-GIM_AllTrig\") ]; + vsX="EK"; % ["Ek"/"En"/"Energy","mm"/"r"/"range","ID"/"IDs"] end %% check of user input data @@ -70,6 +71,7 @@ error("please specify a label for each data set"); end end +if (~exist("vsX","var")), vsX="mm"; end %% clear storage % - clear summary data @@ -141,12 +143,19 @@ end %% show data -addIndex=mmsProf; -addLabel="Range [mm]"; -% addIndex=EksProf; -% addLabel="E_k [MeV/u]"; -% addIndex=repmat((1:(size(profiles,2)-1))',[1 size(profiles,4)]); -% addLabel="ID"; +switch upper(vsX) + case {"EK","EN","ENERGY"} + addIndex=EksProf; + addLabel="E_k [MeV/u]"; + case {"ID","IDS"} + addIndex=repmat((1:(size(profiles,2)-1))',[1 size(profiles,4)]); + addLabel="ID"; + case {"MM","R","RANGE"} + addIndex=mmsProf; + addLabel="Range [mm]"; + otherwise + error("Cannot recognise what you want as X-axis in summary overviews: %s!",vsX); +end if (exist("shifts","var")) for iDataAcq=1:nDataSets addIndex(:,iDataAcq)=addIndex(:,iDataAcq)+shifts(iDataAcq); @@ -168,8 +177,16 @@ CompBars=BARsSumm(:,:,iDataSumm); CompBars(:,:,2)=BARsProf(:,:,iDataAcq); CompFwhms=FWHMsSumm(:,:,iDataSumm); CompFwhms(:,:,2)=FWHMsProf(:,:,iDataAcq); CompInts=INTsSumm(:,:,iDataSumm); CompInts(:,:,2)=INTsProf(:,:,iDataAcq); -% CompXs=mmsSumm(:,iDataAcq); CompXs(:,2)=mmsProf(:,iDataAcq); - CompXs=(1:size(BARsSumm,1))'; CompXs(:,2)=addIndex(:,iDataAcq); + switch upper(vsX) + case {"EK","EN","ENERGY"} + CompXs=EksSumm(:,iDataAcq); CompXs(:,2)=EksProf(:,iDataAcq); + case {"ID","IDS"} + CompXs=(1:size(BARsSumm,1))'; CompXs(:,2)=addIndex(:,iDataAcq); + case {"MM","R","RANGE"} + CompXs=mmsSumm(:,iDataAcq); CompXs(:,2)=mmsProf(:,iDataAcq); + otherwise + error("Cannot recognise what you want as X-axis in summary overviews: %s!",vsX); + end ShowBeamProfilesSummaryData(CompBars,CompFwhms,CompInts,missing(),CompXs,addLabel,... ["summary data" "stat on profiles"],missing(),sprintf("%s - %s - summary vs profile stats",myTit,myLabels(iDataAcq))); end diff --git a/MADX-tracking/ReadStartingConditions.m b/MADX-tracking/StartConditionsRead.m similarity index 100% rename from MADX-tracking/ReadStartingConditions.m rename to MADX-tracking/StartConditionsRead.m From f7c709976f5ff9401c84bd4c87bffea7f766bb08 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 5 Jul 2023 14:15:00 +0200 Subject: [PATCH 05/26] Auto-detection of S-coordinate (centre/end of element) when plotting the lattice --- MADX-optics/PlotLattice.m | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/MADX-optics/PlotLattice.m b/MADX-optics/PlotLattice.m index 659dca4..93dbbe7 100644 --- a/MADX-optics/PlotLattice.m +++ b/MADX-optics/PlotLattice.m @@ -41,7 +41,13 @@ function PlotLattice(geometry) maxS=max(geometry{colS})+geometry{colL}(end); minS=min(geometry{colS}); - + + % auto-detect if S-coordinate is given at the end or at centre of + % element + iL=find(geometry{colL}>0); + iL=iL(1); + lEnd=(geometry{colS}(iL)==geometry{colS}(iL-1)+geometry{colL}(iL)); + % horizontal line, above which elements are displayed plot([minS maxS],[0 0],'k'); @@ -52,7 +58,7 @@ function PlotLattice(geometry) for ii = 1:length(indices) % get position along line and extension - [s,ds]=getSdS(geometry{colS}(indices(ii)),geometry{colL}(indices(ii))); + [s,ds]=getSdS(geometry{colS}(indices(ii)),geometry{colL}(indices(ii)),lEnd); % position of the box if ( length(yPosBoxes{jj}) == 2 ) @@ -79,7 +85,7 @@ function PlotLattice(geometry) % KICKERs: special treatement indices=find(contains(geometry{colKey},'KICKER')); for ii = 1:length(indices) - [s,ds]=getSdS(geometry{colS}(indices(ii)),geometry{colL}(indices(ii))); + [s,ds]=getSdS(geometry{colS}(indices(ii)),geometry{colL}(indices(ii)),lEnd); if contains(geometry{colName}(indices(ii)),'SP') % septum setBox(s,ds,-0.5,1,'m','k'); @@ -106,8 +112,13 @@ function PlotLattice(geometry) end -function [z,dz]=getSdS(S,L) - z=S-L; +function [z,dz]=getSdS(S,L,lEnd) + if (~exist("lEnd","var")), lEnd=true; end % S-coord given at end + if (lEnd) + z=S-L; + else + z=S-L/2.; + end dz=L; end From 4e5e7e21b38fd8ee220f68166775b7dc65a642cb Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 5 Jul 2023 14:16:06 +0200 Subject: [PATCH 06/26] A bit of code refactoring: * `MADX-optics/CompareOptics.m`: using a switch MatLab statement instead of a chain of if-elseif-else statements * `MADX-optics/ParseTfsTableHeader.m`: smarter parsing of .tfs file header * `MADX-tracking/PlotLossMap.m` and `MADX-tracking/ShowLossMap.m`: indices are used only when needed * `MADX-tracking/ReadLosses.m`: smarter handling of optional arguments to function interface + file format determined by `GetColumnsAndMappingTFS` * `MADX-tracking/ScatterPlot.m`: arrangement of plots determined by `GetNrowsNcols`, column names and columns to be plotted determined by `GetColumnsAndMappingTFS`, indices as used only when needed; * `general/Plot2DHistograms.m`: start showing 3D histograms as 2D color maps; --- MADX-optics/CompareOptics.m | 31 +++++++------ MADX-optics/ParseTfsTableHeader.m | 75 +++++++++---------------------- MADX-tracking/PlotLossMap.m | 12 +++-- MADX-tracking/ReadLosses.m | 69 ++++++++++++++-------------- MADX-tracking/ScatterPlot.m | 34 +++++++++----- general/Plot2DHistograms.m | 1 + 6 files changed, 107 insertions(+), 115 deletions(-) diff --git a/MADX-optics/CompareOptics.m b/MADX-optics/CompareOptics.m index c37aa29..05dd98f 100644 --- a/MADX-optics/CompareOptics.m +++ b/MADX-optics/CompareOptics.m @@ -11,20 +11,23 @@ function CompareOptics(optics,labels,geometry,whats,myTitle,emig,sigdpp) uppWhats=upper(whats); switch length(uppWhats) case 1 - if ( strcmp(uppWhats,"BET" ) ) - usrWhats=[ "BETX" "BETY" ]; - elseif ( strcmp(uppWhats,"D" ) ) - usrWhats=[ "DX" "DY" ]; - elseif ( strcmp(uppWhats,"CO" ) ) - usrWhats=[ "X" "Y" ]; - elseif ( strcmp(uppWhats,"SIG" ) ) - usrWhats=[ "SIGX" "SIGY" ]; - elseif ( strcmp(uppWhats,"SIGP" ) ) - usrWhats=[ "SIGPX" "SIGPY" ]; - elseif ( strcmp(uppWhats,"ENV" ) ) - usrWhats=[ "ENVX" "ENVY" ]; - else - error("cannot undestand single what %s!",whats); + switch uppWhats + case "BET" + usrWhats=[ "BETX" "BETY" ]; + case "MU" + usrWhats=[ "MUX" "MUY" ]; + case "D" + usrWhats=[ "DX" "DY" ]; + case "CO" + usrWhats=[ "X" "Y" ]; + case "SIG" + usrWhats=[ "SIGX" "SIGY" ]; + case "SIGP" + usrWhats=[ "SIGPX" "SIGPY" ]; + case "ENV" + usrWhats=[ "ENVX" "ENVY" ]; + otherwise + error("cannot undestand single what %s!",whats); end case 2 usrWhats=uppWhats; diff --git a/MADX-optics/ParseTfsTableHeader.m b/MADX-optics/ParseTfsTableHeader.m index 9ed7477..123dda4 100644 --- a/MADX-optics/ParseTfsTableHeader.m +++ b/MADX-optics/ParseTfsTableHeader.m @@ -1,5 +1,5 @@ function [Qx,Qy,DQx,DQy,Laccel,headerNames,headerValues] = ... - ParseTfsTableHeader(fileNames,nHeader) + ParseTfsTableHeader(fileNames) % ParseTfsTableHeader read header of TFS file. % % The function returns all the names and values of the parameters listed @@ -24,65 +24,34 @@ % headerNames,headerValues: arrays storing field name and value of each % header line. Returned for user convenience; % -% optional input arguments: -% nHeader: number of header lines; -% if not specified, the default value (i.e. 48) is used; -% Nota Bene: the column names and formats are considered as part of the -% header; -% % See also ParseTfsTable. - nHeaderUsr=48; - if ( exist('nHeader','var') ) - nHeaderUsr=nHeader; - end - switch nHeaderUsr - case{48} - % first 4 lines of general text - iStart01=1; - iStop01=4; - % part of header with actual data - iStart02=5; - iStop02=42; - % additional 4 lines of general text - iStart03=43; - iStop03=46; - otherwise - error('unknown number of header lines!'); - return - end - if ( length(fileNames)==1 ) - fid = fopen(fileNames,'r'); + myGets=[ "Q1" "Q2" "DQ1" "DQ2" "LENGTH" ]; + myNums=NaN(size(myGets)); headerNames=[]; headerValues=[]; - for ii = iStart01:iStop01 - % text data - tmp = textscan(fgetl(fid),'@ %s %s %s'); - % headerNames=[headerNames tmp{1}]; - % headerValues=[headerValues tmp{3}]; - end - for ii = iStart02:iStop02 - % real data - tmpLine=fgetl(fid); - tmp = textscan(tmpLine,'@ %s %s %f'); - headerNames=[headerNames tmp{1}]; - headerValues=[headerValues tmp{3}]; - end - for ii = iStart03:iStop03 - % text data - tmpLine=fgetl(fid); - tmp = textscan(tmpLine,'@ %s %s %s'); - % headerNames=[headerNames tmp{1}]; - % headerValues=[headerValues tmp{2}]; + fid = fopen(fileNames,'r'); + while (~feof(fid)) + tmpL=fgetl(fid); + if (startsWith(tmpL,"@")) + tmp = textscan(tmpL,'@ %s %s %s'); + myMatch = strcmpi(myGets,tmp{1}); + if (any(myMatch)), myNums(myMatch)=str2double(tmp{3}); end + if (~startsWith(tmp{3},'"')) + headerNames=[ headerNames tmp{1} ]; + headerValues=[ headerValues tmp{3} ]; + end + else + break + end end fclose(fid); - - Qx=headerValues(find(strcmp(headerNames,'Q1'))); - Qy=headerValues(find(strcmp(headerNames,'Q2'))); - DQx=headerValues(find(strcmp(headerNames,'DQ1'))); - DQy=headerValues(find(strcmp(headerNames,'DQ2'))); - Laccel=headerValues(find(strcmp(headerNames,'LENGTH'))); + Qx=myNums(strcmpi(myGets,"Q1")); + Qy=myNums(strcmpi(myGets,"Q2")); + DQx=myNums(strcmpi(myGets,"DQ1")); + DQy=myNums(strcmpi(myGets,"DQ2")); + Laccel=myNums(strcmpi(myGets,"LENGTH")); else Qx=[]; Qy=[]; diff --git a/MADX-tracking/PlotLossMap.m b/MADX-tracking/PlotLossMap.m index dcc86d9..251e52a 100644 --- a/MADX-tracking/PlotLossMap.m +++ b/MADX-tracking/PlotLossMap.m @@ -1,15 +1,15 @@ -function [hh,cc] = PlotLossMap(losses,indeces,dS,maxS,minS,Nbins) +function [hh,cc] = PlotLossMap(losses,indices,dS,maxS,minS,Nbins) % PlotLossMap plot the longitudinal distribution of losses (histogram) % % Nota bene: the function does not create a figure or a subplot, they % should be created beforehand; % -% [hh,cc] = PlotLossMap(losses,indeces,dS,maxS,minS,Nbins) +% [hh,cc] = PlotLossMap(losses,indices,dS,maxS,minS,Nbins) % % input arguments: % losses: table of losses. For the format of the table, please see % GetVariablesAndMappingParticleData -% indeces: index of particles to be plotted (eg to apply some filtering); +% indices: index of particles to be plotted (eg to apply some filtering); % % optional input arguments: % dS: bin width [m]; @@ -46,7 +46,11 @@ % compute histogram edges=[minSUsr:dsUsr:maxSUsr]; - [hh,cc] = histcounts(losses{colS}(indeces),edges); + if (ismissing(indices)) + [hh,cc] = histcounts(losses{colS},edges); + else + [hh,cc] = histcounts(losses{colS}(indices),edges); + end edges = cc(2:end) - 0.5*(cc(2)-cc(1)); % plot histogram diff --git a/MADX-tracking/ReadLosses.m b/MADX-tracking/ReadLosses.m index 63140f5..4fa6671 100644 --- a/MADX-tracking/ReadLosses.m +++ b/MADX-tracking/ReadLosses.m @@ -1,11 +1,9 @@ -function [ losses, nLosses ] = ReadLosses(path2Files,nInitial,delays,nExpected) -% READLOSSES Read losses from MAD-X tracking. +function [ losses, nTotLosses ] = ReadLosses(path2Files,nInitial,delays,nExpected) +% READLOSSES Read losses from MAD-X/PTC tracking. % % example format of line of parsed files: % NUMBER, TURN, X, PX, Y, PY, T, PT, S, E, ELEMENT % -% [ losses, nLosses ] = ReadLosses(path2Files,nInitial,delays) -% % input arguments: % path2Files: path to files to parse; % example: "delay_%02u\\beam2track*\\losses.tfs" @@ -20,53 +18,56 @@ % % output arguments: % losses: cell array of particle coordinates; -% nLosses: total number of particles read +% nTotLosses: total number of particles read - losses=[]; - nLosses=0; - if ( exist('delays','var') ) - kMax=length(delays); + if (~exist('nInitial','var')), nInitial=missing(); end + if (~exist('delays','var')), delays=missing(); end + if (~exist('nExpected','var')), nExpected=missing(); end + + if (ismissing(nInitial)), cumInitial=0; else cumInitial=cumsum(nInitial); end + if (ismissing(delays)) + kMax=length(path2Files); else - kMax=1; - end - if ( exist('nInitial','var') ) - cumInitial=cumsum(nInitial); + kMax=length(delays); + if (length(path2Files)~=kMax) + error("Inconsistent number of delays (%d) and paths where to find files (%d)!", ... + kMax,length(path2Files)); + end end + + [ colNames, colUnits, colFacts, mapping, readFormat ] = ... + GetColumnsAndMappingTFS("losses"); + + losses=[]; + nTotLosses=0; for kk=1:kMax - if ( exist('delays','var') ) - files=dir(sprintf(path2Files,delays(kk))); - else - files=dir(path2Files); - end + files=dir(path2Files(kk)); for ii=1:length(files) lossesFileName=sprintf("%s\\%s",files(ii).folder,files(ii).name); fprintf("reading file %s ...\n",lossesFileName); fileID = fopen(lossesFileName,'r'); - tmpLosses = textscan(fileID,'%f %f %f %f %f %f %f %f %f %f %s','HeaderLines',8); + tmpLosses = textscan(fileID,readFormat,'HeaderLines',8); fclose(fileID); - if ( exist('delays','var') ) - tmpLosses{2}=tmpLosses{2}+delays(kk); + if (~ismissing(delays)) + iCol=mapping(strcmpi(colNames,"TURN")); + tmpLosses{iCol}=tmpLosses{iCol}+delays(kk); end - if ( exist('nInitial','var') ) - if ( ii>1 ) - tmpLosses{1}=tmpLosses{1}+cumInitial(ii-1); - end + if (~ismissing(nInitial) && ii>1) + iCol=mapping(strcmpi(colNames,"NUMBER")); + tmpLosses{iCol}=tmpLosses{iCol}+cumInitial(ii-1); end - nLosses=nLosses+length(tmpLosses{1}); + nLosses=length(tmpLosses{1}); if (length(losses)==0) losses=tmpLosses; else - for jj=1:length(losses) - losses{jj}=[ losses{jj} ; tmpLosses{jj} ]; - end + losses{nTotLosses+1:nTotLosses+nLosses}=tmpLosses; end + nTotLosses=nTotLosses+nLosses; end end - fprintf("...acquired %u particle losses; \n",nLosses); + fprintf("...acquired %u particle losses; \n",nTotLosses); - if ( exist('nExpected','var') ) - if ( nLosses > nExpected ) - warning('...something wrong: I was expecting at most %u points whereas I have parsed %u',nExpected,nLosses); - end + if ( ~ismissing(nExpected) && nTotLosses>nExpected ) + warning('...something wrong: I was expecting at most %u points whereas I have parsed %u',nExpected,nTotLosses); end end \ No newline at end of file diff --git a/MADX-tracking/ScatterPlot.m b/MADX-tracking/ScatterPlot.m index 8274262..69c1d9f 100644 --- a/MADX-tracking/ScatterPlot.m +++ b/MADX-tracking/ScatterPlot.m @@ -1,13 +1,13 @@ -function ScatterPlot(data,indeces,xIDs,yIDs,whichData,currTitle) +function ScatterPlot(data,indices,xIDs,yIDs,whichData,currTitle) % ScatterPlot Prepare scatter plots % % The user decides what to show % -% ScatterPlot(data,indeces,xIDs,yIDs,whichData,currTitle) +% ScatterPlot(data,indices,xIDs,yIDs,whichData,currTitle) % % input arguments: % data: table with particle coordinates; -% indeces: index of particles to be plotted (eg to apply some filtering); +% indices: index of particles to be plotted (eg to apply some filtering); % xIDs, yIDs: ID of the quantity to be plotted: % whichData: format of the table; % currTitle: title of the plot; @@ -18,17 +18,31 @@ function ScatterPlot(data,indeces,xIDs,yIDs,whichData,currTitle) error('...number of x-axes different from number of y-axes!'); return end - nn=ceil(sqrt(length(xIDs))); + [nRows,nCols]=GetNrowsNcols(length(xIDs)); - [ colNames, colUnits, colFacts, mapping ] = GetVariablesAndMapping(whichData); + [ colNames, colUnits, colFacts, mapping ] = GetColumnsAndMappingTFS(whichData); actualTitle=sprintf('%s - Scatter plot - %s',whichData,currTitle); ff=figure('Name',actualTitle,'NumberTitle','off'); - for ii=1:size(xIDs,2) - axt=subplot(nn,nn,ii); - plot(data{mapping(xIDs(ii))}(indeces)*colFacts(xIDs(ii)),data{mapping(yIDs(ii))}(indeces)*colFacts(yIDs(ii)),'r*'); - xlabel(sprintf('%s [%s]',colNames(xIDs(ii)),colUnits(xIDs(ii)))); - ylabel(sprintf('%s [%s]',colNames(yIDs(ii)),colUnits(yIDs(ii)))); + for ii=1:length(xIDs) + iColX=mapping(strcmpi(colNames,xIDs(ii))); + iColY=mapping(strcmpi(colNames,yIDs(ii))); + axt=subplot(nRows,nCols,ii); + if (ismissing(indices)) + if (iscell(data)) + plot(data{iColX}*colFacts(iColX),data{iColY}*colFacts(iColY),'r.'); + else + plot(data(:,iColX)*colFacts(iColX),data(:,iColY)*colFacts(iColY),'r.'); + end + else + if (iscell(data)) + plot(data{iColX}(indices)*colFacts(iColX),data{iColY}(indices)*colFacts(iColY),'r.'); + else + plot(data(indices,iColX)*colFacts(iColX),data(indices,iColY)*colFacts(iColY),'r.'); + end + end + xlabel(sprintf('%s [%s]',colNames(iColX),colUnits(iColX))); + ylabel(sprintf('%s [%s]',colNames(iColY),colUnits(iColY))); grid on; end sgtitle(actualTitle); diff --git a/general/Plot2DHistograms.m b/general/Plot2DHistograms.m index 013f9fd..0fdb4bb 100644 --- a/general/Plot2DHistograms.m +++ b/general/Plot2DHistograms.m @@ -7,6 +7,7 @@ if ( lHist ) hh=histogram2('XBinEdges',xShow,'YBinEdges',yShow,... 'BinCounts',showMe,'DisplayStyle','tile','ShowEmptyBins','on'); + view(2); % starts appearing as a 2D plot else plot(showMe(:,1),showMe(:,2),"k."); end From edc0fd283b78476b5839efca70525e208c515a8c Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 5 Jul 2023 14:17:20 +0200 Subject: [PATCH 07/26] Read/Write starting conditions for MADX/PTC tracking --- MADX-tracking/StartConditionsFormat.m | 13 +++++++++++++ MADX-tracking/StartConditionsRead.m | 19 +++++++++++-------- MADX-tracking/StartConditionsWrite.m | 11 +++++++++++ 3 files changed, 35 insertions(+), 8 deletions(-) create mode 100644 MADX-tracking/StartConditionsFormat.m create mode 100644 MADX-tracking/StartConditionsWrite.m diff --git a/MADX-tracking/StartConditionsFormat.m b/MADX-tracking/StartConditionsFormat.m new file mode 100644 index 0000000..b07f8d7 --- /dev/null +++ b/MADX-tracking/StartConditionsFormat.m @@ -0,0 +1,13 @@ +function myFmt=StartConditionsFormat(tracker) +% example format of line of parsed files (argument of ID is read as well): +% START, X=-0.01, PX=-0.01, Y=0, PY=0, T=0, PT=0; ! ID=1 +% PTC_START, X=-0.01, PX=-0.01, Y=0, PY=0, T=0, PT=0; ! ID=1 + switch upper(tracker) + case {"MADX","MAD","MAD-X"} + myFmt='START, X=%f, PX=%f, Y=%f, PY=%f, T=%f, PT=%f;%*[^\n]'; + case "PTC" + myFmt='PTC_START, X=%f, PX=%f, Y=%f, PY=%f, T=%f, PT=%f;%*[^\n]'; + otherwise + error("Unknown tracker %s!",tracker); + end +end diff --git a/MADX-tracking/StartConditionsRead.m b/MADX-tracking/StartConditionsRead.m index a93a39d..a2225d4 100644 --- a/MADX-tracking/StartConditionsRead.m +++ b/MADX-tracking/StartConditionsRead.m @@ -1,14 +1,13 @@ -function [ particles, nInitial, nTot ] = ReadStartingConditions(path2Files,nExpected) -% READSTARTINGCONDITIONS Read starting conditions of MAD-X tracking. +function [ particles, nInitial, nTot ] = StartConditionsRead(path2Files,nExpected,tracker) +% STARTCONDITIONSREAD Read starting conditions of MAD-X/PTC tracking. % % example format of line of parsed files (argument of ID is read as well): -% START, X=-0.01, PX=-0.01, Y=0, PY=0, T=0, PT=0; ! ID=1 +% START, X=-0.01, PX=-0.01, Y=0, PY=0, T=0, PT=0; +% PTC_START, X=-0.01, PX=-0.01, Y=0, PY=0, T=0, PT=0; % for the time being: % no empty lines; % no MAD-X comments; -% a line heading with 'reutrn' is skipped; -% -% [ particles, nInitial, nTot ] = ReadStartingConditions(path2Files,nExpected) +% a line heading with 'return' is skipped; % % input arguments: % path2Files: path to files to parse; @@ -16,12 +15,16 @@ % nExpected: expected number of initial conditions; % if given, a check on the number of conditions parsed is % performed; +% tracker: MADX, PTC % % output arguments: % particles: cell array of particle coordinates; % nInitial: array of cumulative number of particles; % nTot: total number of particles read + if ( ~exist('nExpected','var') ), nExpected=missing(); end + if ( ~exist('tracker','var') ), tracker="MADX"; end + myFmt=StartConditionsFormat(tracker); particles=[]; files=dir(path2Files); nInitial=[]; @@ -29,7 +32,7 @@ tmpFileName=sprintf("%s\\%s",files(ii).folder,files(ii).name); fprintf("reading file %s ...\n",tmpFileName); fileID = fopen(tmpFileName,'r'); - tmpParticles = textscan(fileID,'START, X=%f, PX=%f, Y=%f, PY=%f, T=%f, PT=%f; ! ID=%f','TreatAsEmpty','return;'); + tmpParticles = textscan(fileID,myFmt,'TreatAsEmpty','return;'); nInitial=[ nInitial length(tmpParticles{1}) ]; fclose(fileID); if (length(particles)==0) @@ -42,7 +45,7 @@ end nTot=sum(nInitial); fprintf("...acquired %u starting conditions; \n",nTot); - if ( exist('nExpected','var') ) + if ( ~ismissing(nExpected) ) if ( nTot ~= nExpected ) warning('...something wrong: I was expecting %u points whereas I have parsed %u',nExpected,nTot); else diff --git a/MADX-tracking/StartConditionsWrite.m b/MADX-tracking/StartConditionsWrite.m new file mode 100644 index 0000000..54bce51 --- /dev/null +++ b/MADX-tracking/StartConditionsWrite.m @@ -0,0 +1,11 @@ +function StartConditionsWrite(FileName,particles,tracker) + if ( ~exist('tracker','var') ), tracker="MADX"; end + fprintf("writing starting conditions of %d particles to file %s (%s format)...\n",length(particles),FileName,tracker); + myFmt=StartConditionsFormat(tracker); + fileID = fopen(FileName,'w'); + for jj=1:length(particles) + fprintf(fileID,"%s\n",sprintf(myFmt,particles(jj,:))); + end + fclose(fileID); + fprintf("...done;\n"); +end From 57cf6e602245bb5320729622c14a1d7285c7cb7e Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 5 Jul 2023 14:17:55 +0200 Subject: [PATCH 08/26] Implementing parsing of track/loss files (MADX/PTC) --- MADX-optics/GetColumnsAndMappingTFS.m | 30 ++++++++- MADX-tracking/ReadSurvivors.m | 88 +++++++++++++++++++++++++++ 2 files changed, 115 insertions(+), 3 deletions(-) create mode 100644 MADX-tracking/ReadSurvivors.m diff --git a/MADX-optics/GetColumnsAndMappingTFS.m b/MADX-optics/GetColumnsAndMappingTFS.m index 7a913b2..90a91a5 100644 --- a/MADX-optics/GetColumnsAndMappingTFS.m +++ b/MADX-optics/GetColumnsAndMappingTFS.m @@ -12,6 +12,8 @@ % 'GEOMETRY': lattice geometry by MAD-X; % 'RMATRIX': (selected elements of) response matrix by MAD-X; % 'CURR': a TFS table listing scan currents; +% 'LOSSES': an ASCII table listing particle losses; +% 'TRACKS': an ASCII table listing particle coordinates; % genBy: format of the table: % 'TWISS': a TFS table as generated by a MADX TWISS command; % 'SCAN': a TFS table as generated by a current scan; @@ -38,6 +40,8 @@ % . response matrix: % NAME, KEYWORD, L, S, RE11, RE12, RE21, RE22, RE16, RE26, RE33, RE34, % RE43, RE44, RE36, RE46, RE51, RE52, RE55, RE56, RE66; +% . losses: +% NUMBER, TURN, X, PX, Y, PY, T, PT, S, E, ELEMENT % - generated by a scan (genBy=="SCAN") % . optics: % Brho[Tm], BP[mm], myID[], BETX, ALFX, BETY, ALFY, X, PX, Y, PY, @@ -49,6 +53,8 @@ % Brho[Tm], BP[mm], myID[], RE11, RE12, RE21, RE22, RE16, RE26, RE33, RE34, % RE43, RE44, RE36, RE46, RE51, RE52, RE55, RE56, RE66; % +% NB: current .tfs tables are custom-made, hence no pre-defined format. +% % See also ParseTfsTable. % default genBy if (~exist("genBy","var") | ismissing(genBy)) @@ -92,7 +98,7 @@ error("NO pre-defined format for CURR .tfs files!"); end switch upper(genBy) - case "TWISS" + case {"TWISS","PTC_TWISS"} switch upper(whichTFS) case {'OPT','OPTICS'} colNames=[ "NAME" "KEYWORD" "L" "S" "BETX" "ALFX" "BETY" "ALFY" ... @@ -102,7 +108,7 @@ "m" "" "m" "" "m" "" "m" "" ... "2\pi" "2\pi"]; readFormat = '%s %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f'; - case {'GEO','GEOMETRY'} + case {'GEO','GEOMETRY','GEOM'} colNames=[ "NAME" "KEYWORD" "L" "S" ... "KICK" "HKICK" "VKICK" "ANGLE" "K0L" "K1L" "K2L" ... "APERTYPE" "APER_1" "APER_2" "APER_3" "APER_4" "APOFF_1" "APOFF_2" ]; @@ -110,7 +116,7 @@ "rad" "rad" "rad" "rad" "rad" "rad m^{-1}" "rad m^{-2}" ... "" "m" "m" "m" "m" "m" "m" ]; readFormat = '%s %s %f %f %f %f %f %f %f %f %f %s %f %f %f %f %f %f %f'; - case {'RM','RMATRIX'} + case {'RM','RMATRIX','RE'} colNames=[ "NAME" "KEYWORD" "L" "S" ... "RE11" "RE12" "RE21" "RE22" "RE16" "RE26" ... "RE33" "RE34" "RE43" "RE44" "RE36" "RE46" ... @@ -120,6 +126,22 @@ "" "m rad^{-1}" "rad m^{-1}" "" "m" "rad" ... "s m^{-1}" "s rad^{-1}" "" "s" "m" ]; readFormat = '%s %s %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f %f'; + case {'LOSS','LOSSES','LOST'} + colNames=[ "NUMBER" "TURN" ... + "X" "PX" "Y" "PY" "T" "PT" ... + "S" "E" "ELEMENT" ]; + colUnits=[ "" "" ... + "m" "" "m" "" "m" "" ... + "m" "GeV" "" ]; + readFormat = '%f %f %f %f %f %f %f %f %f %f %s'; + case {'TRACK','TRACKS','TRACKING','SURV','SURVS','SURVIVOR','SURVIVORS'} + colNames=[ "NUMBER" "TURN" ... + "X" "PX" "Y" "PY" "T" "PT" ... + "S" "E" ]; + colUnits=[ "" "" ... + "m" "" "m" "" "m" "" ... + "m" "GeV" ]; + readFormat = '%f %f %f %f %f %f %f %f %f %f'; otherwise error('which column mapping for TFS table?'); end @@ -165,6 +187,8 @@ '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep,'%f' ); + case {'LOSS','LOSSES','TRACK','TRACKS','TRACKING'} + error('... %s .tfs type generated by %s NOT available yet!',whichTFS,genBy); otherwise error('which column mapping for TFS table?'); end diff --git a/MADX-tracking/ReadSurvivors.m b/MADX-tracking/ReadSurvivors.m new file mode 100644 index 0000000..d1aabe1 --- /dev/null +++ b/MADX-tracking/ReadSurvivors.m @@ -0,0 +1,88 @@ +function [ particles, obsPoints, nParts ] = ReadSurvivors(path2Files) +% READSURVIVORS Read survivors from MAD-X/PTC tracking. +% +% example format of line of parsed files (argument of ID is read as well): +% * NUMBER TURN X PX Y PY T PT S E +% #segment 1 32 200000 0 start +% 1 0 -0.02 -0.002 -0.013 -0.006 0 -0.0013 0 0 +% +% input arguments: +% path2Files: path to files to parse; +% +% output arguments: +% particles: table of particle coordinates; +% the third dimension identifies the obervation point; +% obsPoints: list of names of the obervation points; + + [ colNames, colUnits, colFacts, mapping, readFormat ] = ... + GetColumnsAndMappingTFS("track"); + + % acquire files + particles=missing(); obsPoints=missing(); nParts=missing(); + for kk=1:length(path2Files) + files=dir(path2Files(kk)); + for ii=1:length(files) + tmpFileName=sprintf("%s\\%s",files(ii).folder,files(ii).name); + fprintf("reading file %s ...\n",tmpFileName); + fileID = fopen(tmpFileName,'r'); + CO6D=zeros(6,1); CO6Dnames=["XC" "PXC" "YC" "PYC" "TC" "PTC"]; + while (~feof(fileID)) + tLine=strip(fgetl(fileID)); + if (startsWith(tLine,"@")) + % header + myFields=split(tLine); + iCol=strcmpi(CO6Dnames,myFields{2}); + if (any(iCol)) + CO6D(iCol)=str2double(myFields{2}); + fprintf("...found %s=%f in header,\n",CO6Dnames(iCol),CO6D(iCol)); + end + elseif (startsWith(tLine,"#")) + % new section + myFields=split(tLine); + nParts=ExpandMat(nParts,double(string(myFields{4}))); + clear tmpParticles; + tmpParticles=cell2mat(textscan(fileID,readFormat)); + fprintf("...acquired %d particles at %s;\n",nParts(end),myFields{end}); + % add closed orbit + if (strcmpi(myFields{end},"START")) + for jj=1:length(CO6Dnames) + if (CO6D(jj)~=0.0) + iCol=mapping(strcmpi(colNames,extractBetween(CO6Dnames(ii),1,strlength(CO6Dnames(ii))-1))); + tmpParticles(:,iCol)=tmpParticles(:,iCol)+CO6D(jj); + fprintf("...taking into account 6D closed orbit on %s at START;\n",colNames(iCol)); + end + end + end + % store in memory + particles=ExpandMat(particles,tmpParticles); + obsPoints=ExpandMat(obsPoints,string(myFields{end})); + end + end + fclose(fileID); + end + end + fprintf("...done with parsing files;\n"); + + % compact storage + uObsPoints=unique(obsPoints); + if (length(uObsPoints)~=length(obsPoints)) + fprintf("...compacting from %d data sets to %d (unique obs points);\n",length(obsPoints),length(uObsPoints)); + readParticles=particles; readObsPoints=obsPoints; readNparts=nParts; + clear particles obsPoints nParts; + particles=missing(); obsPoints=missing(); nParts=missing(); + for ii=1:length(uObsPoints) + iObs=find(strcmp(readObsPoints,uObsPoints(ii))); + clear tmpParticles; + nParts=ExpandMat(nParts,sum(readNparts(iObs))); + tmpParticles=NaN(nParts(end),length(colNames)); + jRead=0; + for jj=1:length(iObs) + tmpParticles(jRead+1:jRead+readNparts(iObs(jj)),:)=readParticles(1:readNparts(iObs(jj)),:,iObs(jj)); + jRead=jRead+readNparts(iObs(jj)); + end + particles=ExpandMat(particles,tmpParticles); + obsPoints=ExpandMat(obsPoints,uObsPoints(ii)); + end + fprintf("...done with compacting data storage;\n"); + end +end \ No newline at end of file From 12a3b15b1618b99f40a257d4f2bd67b0a7e5466e Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 5 Jul 2023 14:18:58 +0200 Subject: [PATCH 09/26] forgotten file in a previous code-refactoring commit --- MADX-tracking/ShowLossMap.m | 24 ++++++++++++++++-------- 1 file changed, 16 insertions(+), 8 deletions(-) diff --git a/MADX-tracking/ShowLossMap.m b/MADX-tracking/ShowLossMap.m index 136f165..f24b0ac 100644 --- a/MADX-tracking/ShowLossMap.m +++ b/MADX-tracking/ShowLossMap.m @@ -1,14 +1,14 @@ -function [hh,cc]=ShowLossMap(losses,indeces,dS,geometry,allAperH,allAperV) +function [hh,cc]=ShowLossMap(losses,indices,dS,geometry,allAperH,allAperV) % ShowLossMap embed the longitudinal distribution of losses (histogram) % into a nice plot, possibly with lattice structure and % aperture model % -% [hh,cc] = ShowLossMap(losses,indeces,dS,geometry,allAperH,allAperV) +% [hh,cc] = ShowLossMap(losses,indices,dS,geometry,allAperH,allAperV) % % input arguments: % losses: table of losses. For the format of the table, please see % GetVariablesAndMappingParticleData -% indeces: index of particles to be plotted (eg to apply some filtering); +% indices: index of particles to be plotted (eg to apply some filtering); % % optional input arguments: % dS: bin width [m]; @@ -73,8 +73,12 @@ PlotAperture(allAperH{1},allAperH{2},allAperH{3},allAperH{4},allAperH{5}); hold on; % losses - plot(losses{colS}(indeces),losses{colX}(indeces)*1000,'r*'); - title(sprintf('H plane')); + if (ismissing(indices)) + plot(losses{colS},losses{colX},'r.'); + else + plot(losses{colS}(indices),losses{colX}(indices),'r.'); + end + title(sprintf('H plane')); grid on; xlim([minS maxS]); end @@ -86,8 +90,12 @@ PlotAperture(allAperV{1},allAperV{2},allAperV{3},allAperV{4},allAperV{5}); hold on; % losses - plot(losses{colS}(indeces),losses{colY}(indeces)*1000,'r*'); - title(sprintf('V plane')); + if (ismissing(indices)) + plot(losses{colS},losses{colY},'r.'); + else + plot(losses{colS}(indices),losses{colY}(indices),'r.'); + end + title(sprintf('V plane')); grid on; xlim([minS maxS]); end @@ -95,7 +103,7 @@ iPlot=iPlot+1; tmpAx=subplot(nPlots,1,iPlot); axs=[ axs tmpAx ]; - [hh,cc] = PlotLossMap(losses,indeces,dsUsr,maxS,minS); + [hh,cc] = PlotLossMap(losses,indices,dsUsr,maxS,minS); set(tmpAx, 'YScale', 'log'); % manual treatment of yticklabels... % set(tmpAx, 'YTick', [min(ylim):10:max(ylim)]); From 3dc051ddb159e93618fd1675d2f2a21eb84ffd36 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Fri, 7 Jul 2023 16:52:53 +0200 Subject: [PATCH 10/26] Improving the suite of tools for generating/handling simulated beams --- general/Plot2DHistograms.m | 41 ++++++++----- optics/BeamContour.m | 55 ++++++++++++++++++ optics/BeamContourStats.m | 26 +++++++++ optics/BeamNormalise.m | 29 ++++++++++ optics/BeamSample.m | 57 ++++++++++++++++++ .../{BeamSample_Gauss.m => BeamSampleGauss.m} | 2 +- .../{BeamSample_Rect.m => BeamSampleRect.m} | 2 +- optics/BeamStats.m | 34 +++++++++++ optics/BeamStatsShow.m | 58 +++++++++++++++++++ 9 files changed, 286 insertions(+), 18 deletions(-) create mode 100644 optics/BeamContour.m create mode 100644 optics/BeamContourStats.m create mode 100644 optics/BeamNormalise.m create mode 100644 optics/BeamSample.m rename optics/{BeamSample_Gauss.m => BeamSampleGauss.m} (96%) rename optics/{BeamSample_Rect.m => BeamSampleRect.m} (94%) create mode 100644 optics/BeamStats.m create mode 100644 optics/BeamStatsShow.m diff --git a/general/Plot2DHistograms.m b/general/Plot2DHistograms.m index 0fdb4bb..588c3fd 100644 --- a/general/Plot2DHistograms.m +++ b/general/Plot2DHistograms.m @@ -1,6 +1,21 @@ -function [axM,axX,axY]=Plot2DHistograms(showMe,showMe1DX,showMe1DY,xShow,yShow,xShowLabel,yShowLabel,contours,lHist,lSquared) +function [axM,axX,axY]=Plot2DHistograms(showMe,showMe1DX,showMe1DY,xINs,yINs,xShowLabel,yShowLabel,contours,lHist,lSquared,lBinEdges) if ( ~exist("lHist","var") ), lHist=true; end if ( ~exist("lSquared","var") ), lSquared=false; end + if ( ~exist("lBinEdges","var") ), lBinEdges=true; end + + contourFmts=["o" "*"]; + + if (lBinEdges) + xShow=xINs; + yShow=yINs; + xCentres=xINs(2:end)-0.5*diff(xINs(1:2)); + yCentres=yINs(2:end)-0.5*diff(yINs(1:2)); + else + xShow=[xINs-0.5*diff(xINs(1:2)) xINs(end)+0.5*diff(xINs(1:2))]; + yShow=[yINs-0.5*diff(yINs(1:2)) yINs(end)+0.5*diff(yINs(1:2))]; + xCentres=xINs; + yCentres=yINs; + end %% 2D histogram axM=subplot('Position', [0.10, 0.10, 0.6, 0.6]); @@ -14,7 +29,9 @@ if ( exist("contours","var") ) if ( ~all(ismissing(contours),"all") ) for ii=1:size(contours,3) - hold on; plot(contours(:,1,ii),contours(:,2,ii),".-"); + hold on; contourFmt=".-"; + if (ii<=length(contourFmts)), contourFmt=contourFmts(ii); end + plot(contours(:,1,ii),contours(:,2,ii),strcat(contourFmt,"-")); end end end @@ -39,25 +56,17 @@ %% hor variable axX=subplot('Position', [0.10, 0.75, 0.6, 0.15]); - edges = xShow(2:end) - 0.5*(xShow(2)-xShow(1)); - bar(edges,showMe1DX,1); - if ( ~lSquared ) - xlim([min(xShow) max(xShow)]); - end - ylabel('[]'); + bar(xCentres,showMe1DX,1,"EdgeColor","none"); + if ( ~lSquared ), xlim([min(xShow) max(xShow)]); end + ylabel('[]'); grid on; set(axX,'xticklabel',{[]}) - grid on; %% ver variable axY=subplot('Position', [0.75, 0.10, 0.15, 0.6]); - edges = yShow(2:end) - 0.5*(yShow(2)-yShow(1)); - barh(edges,showMe1DY,1); - if ( ~lSquared ) - ylim([min(yShow) max(yShow)]); - end - xlabel('[]'); + barh(yCentres,showMe1DY,1,"EdgeColor","none"); + if ( ~lSquared ), ylim([min(yShow) max(yShow)]); end + xlabel('[]'); grid on; set(axY,'yticklabel',{[]}) - grid on; %% general stuff linkaxes([axM,axX],'x'); diff --git a/optics/BeamContour.m b/optics/BeamContour.m new file mode 100644 index 0000000..a06d308 --- /dev/null +++ b/optics/BeamContour.m @@ -0,0 +1,55 @@ +function [Contours]=BeamContour(dType,alfa,beta,emiGeo,iNorm,bb,hh) +% +% Contours=float(nPoints,2,nPlanes); +% NB: iNorm,bb,hh only for bar-of-charge (optional arguments) +% + fprintf("Contouring beam...\n"); + + %% set up + % - general beam sampling + nPlanes=length(dType); + if (nPlanes~=length(alfa)) + error("Incosistent number of types (%d) and alfa-optics functions (%d) (including NaNs)!",nPlanes,length(alfa)); + end + if (nPlanes~=length(beta)) + error("Incosistent number of types (%d) and beta-optics functions (%d) (including NaNs)!",nPlanes,length(beta)); + end + if (nPlanes~=length(emiGeo)) + error("Incosistent number of types (%d) and geometric emittances (%d) (including NaNs)!",nPlanes,length(emiGeo)); + end + % - specific to BAR-of-charge + nBars=sum(strcmpi(dType,"BAR")); + if (~exist("iNorm","var")), iNorm=false(nBars,1); end + if (~exist("bb","var")), bb=NaN(nBars,1); end + if (~exist("hh","var")), hh=NaN(nBars,1); end + if (nBars~=length(iNorm)) + error("Incosistent number of BARs (%d) and specs for norm/phys units (%d)!",nBars,length(iNorm)); + end + if (nBars~=length(bb)) + error("Incosistent number of BARs (%d) and bar-of-charge base values (%d)!",nBars,length(bb)); + end + if (nBars~=length(hh)) + error("Incosistent number of BARs (%d) and bar-of-charge height values (%d)!",nBars,length(hh)); + end + + %% actually contour + Contours=missing(); + for iPlane=1:nPlanes + switch upper(dType(iPlane)) + case {"B","BAR","BAR-OF-CHARGE"} + Contours=ExpandMat(Contours,GenPointsAlongRectangle(bb(iBar),hh(iBar))); + if (iNorm(iBar)) + Contours(:,:,iPlane)=Norm2Phys(Contours(:,:,iPlane),beta(iPlane),alfa(iPlane),emiGeo(iPlane)); + end + case {"D","DELTA","P","POINT"} + Contours=ExpandMat(Contours,zeros(1,2)); + case {"G","GAUSS","GAUSSIAN"} + Contours=ExpandMat(Contours,GenPointsAlongEllypse(alfa(iPlane),beta(iPlane),emiGeo(iPlane))); + otherwise + error("Unknown distribution type %s!",dType(iPlane)); + end + end + + %% done + fprintf("...done;\n"); +end diff --git a/optics/BeamContourStats.m b/optics/BeamContourStats.m new file mode 100644 index 0000000..df3afa9 --- /dev/null +++ b/optics/BeamContourStats.m @@ -0,0 +1,26 @@ +function [Contours]=BeamContourStats(particles,emiMax) +% particles=float(nParticles,2,nPlanes); + fprintf("Finding fit ellypse on beam...\n"); + + %% set up + nPlanes=size(particles,3); + if (~exist("emiMax","var")), emiMax=false; end + if (emiMax) + fprintf("...based on largest single-particle geometric emittance!\n"); + else + fprintf("...based on RMS single-particle geometric emittance!\n"); + end + + %% actually contour + Contours=missing(); + for iPlane=1:nPlanes + [alpha,beta,emiG]=GetOpticsFromSigmaMatrix(particles(:,:,iPlane)); % ellypse orientation + if (emiMax) + emiG=max(GetSinglePartEmittance(particles(:,:,iPlane),alpha,beta)); % max + end + Contours=ExpandMat(Contours,GenPointsAlongEllypse(alpha,beta,emiG)); + end + + %% done + fprintf("...done;\n"); +end diff --git a/optics/BeamNormalise.m b/optics/BeamNormalise.m new file mode 100644 index 0000000..e73abc2 --- /dev/null +++ b/optics/BeamNormalise.m @@ -0,0 +1,29 @@ +function [NormCoords]=BeamNormalise(PhysCoords,alfa,beta,emiGeo) +% PhysCoords=float(nParticles,2,nPlanes); + fprintf("Normalising beam...\n"); + + %% set up + nPlanes=size(PhysCoords,3); + % defaults: alfa=0; beta=10 m; emiGeo=1E-6 m rad; + if (~exist("alfa","var")), alfa=zeros(nPlanes,1); end + if (~exist("beta","var")), beta=ones(nPlanes,1)*10; end + if (~exist("emiGeo","var")), emiGeo=ones(nPlanes,1)*1E-6; end + if (nPlanes~=length(alfa)) + error("Incosistent number of planes (%d) and alfa-optics functions (%d) (including NaNs)!",nPlanes,length(alfa)); + end + if (nPlanes~=length(beta)) + error("Incosistent number of planes (%d) and beta-optics functions (%d) (including NaNs)!",nPlanes,length(beta)); + end + if (nPlanes~=length(emiGeo)) + error("Incosistent number of planes (%d) and geometric emittances (%d) (including NaNs)!",nPlanes,length(emiGeo)); + end + + %% actually normalise coordinates + NormCoords=missing(); + for iPlane=1:nPlanes + NormCoords=ExpandMat(NormCoords,Phys2Norm(PhysCoords(:,:,iPlane),beta(iPlane),alfa(iPlane),emiGeo(iPlane))); + end + + %% done + fprintf("...done;\n"); +end diff --git a/optics/BeamSample.m b/optics/BeamSample.m new file mode 100644 index 0000000..566adc8 --- /dev/null +++ b/optics/BeamSample.m @@ -0,0 +1,57 @@ +function [beamParts]=BeamSample(dType,nParticles,alfa,beta,emiGeo,iNorm,bb,hh) +% +% beamParts=float(nParticles,2,nPlanes); +% particle coordinates expressed in physical units, i.e. [m] and [rad]! +% NB: iNorm,bb,hh only for bar-of-charge (optional arguments) +% + fprintf("Generating beam of %d particles...\n",nParticles); + + %% set up + % - general beam sampling + nPlanes=length(dType); + if (nPlanes~=length(alfa)) + error("Incosistent number of types (%d) and alfa-optics functions (%d) (including NaNs)!",nPlanes,length(alfa)); + end + if (nPlanes~=length(beta)) + error("Incosistent number of types (%d) and beta-optics functions (%d) (including NaNs)!",nPlanes,length(beta)); + end + if (nPlanes~=length(emiGeo)) + error("Incosistent number of types (%d) and geometric emittances (%d) (including NaNs)!",nPlanes,length(emiGeo)); + end + % - specific to BAR-of-charge + nBars=sum(strcmpi(dType,"BAR")); + if (~exist("iNorm","var")), iNorm=false(nBars,1); end + if (~exist("bb","var")), bb=NaN(nBars,1); end + if (~exist("hh","var")), hh=NaN(nBars,1); end + if (nBars~=length(iNorm)) + error("Incosistent number of BARs (%d) and specs for norm/phys units (%d)!",nBars,length(iNorm)); + end + if (nBars~=length(bb)) + error("Incosistent number of BARs (%d) and bar-of-charge base values (%d)!",nBars,length(bb)); + end + if (nBars~=length(hh)) + error("Incosistent number of BARs (%d) and bar-of-charge height values (%d)!",nBars,length(hh)); + end + + %% actually sample + beamParts=missing(); iBar=0; + for iPlane=1:nPlanes + switch upper(dType(iPlane)) + case {"B","BAR","BAR-OF-CHARGE"} + iBar=iBar+1; + beamParts=ExpandMat(beamParts,BeamSampleRect(nParticles,bb(iBar),hh(iBar))); + if (iNorm(iBar)) + beamParts(:,:,iPlane)=Norm2Phys(beamParts(:,:,iPlane),beta(iPlane),alfa(iPlane),emiGeo(iPlane)); + end + case {"D","DELTA","P","POINT"} + beamParts=ExpandMat(beamParts,zeros(nParticles,2)); + case {"G","GAUSS","GAUSSIAN"} + beamParts=ExpandMat(beamParts,BeamSampleGauss(nParticles,alfa(iPlane),beta(iPlane),emiGeo(iPlane))); + otherwise + error("Unknown distribution type %s!",dType(iPlane)); + end + end + + %% done + fprintf("...done;\n"); +end diff --git a/optics/BeamSample_Gauss.m b/optics/BeamSampleGauss.m similarity index 96% rename from optics/BeamSample_Gauss.m rename to optics/BeamSampleGauss.m index 0d7ed86..0c2df9c 100644 --- a/optics/BeamSample_Gauss.m +++ b/optics/BeamSampleGauss.m @@ -1,4 +1,4 @@ -function MyPoints=BeamSample_Gauss(nPoints,alpha,beta,emiG) +function MyPoints=BeamSampleGauss(nPoints,alpha,beta,emiG) % BeamSample_Gauss(nPoints,alpha,beta,emiG) generates points with a 2D % Gaussian distribution filling the % ellypse described by the given optics diff --git a/optics/BeamSample_Rect.m b/optics/BeamSampleRect.m similarity index 94% rename from optics/BeamSample_Rect.m rename to optics/BeamSampleRect.m index c466426..6b9b676 100644 --- a/optics/BeamSample_Rect.m +++ b/optics/BeamSampleRect.m @@ -1,4 +1,4 @@ -function MyPoints=BeamSample_Rect(nPoints,bb,hh) +function MyPoints=BeamSampleRect(nPoints,bb,hh) % BeamSample_Rect(nPoints,bb,hh) generates points inside an optics % rectangle (e.g. bar of charge), % uniformly distributed; diff --git a/optics/BeamStats.m b/optics/BeamStats.m new file mode 100644 index 0000000..7222178 --- /dev/null +++ b/optics/BeamStats.m @@ -0,0 +1,34 @@ +function [Counts2D,CountsX,CountsY,xb,yb]=BeamStats(particles,what2Analise,declared) +% particles=float(nParticles,2,nPlanes); + fprintf("Beam statistics...\n"); + + %% set up + if (~exist("declared","var")), declared=["X" "PX" "Y" "PY" "T" "PT"]; end + if (ndims(particles)==3) + myParts=reshape(particles,[size(particles,1) size(particles,2)*size(particles,3)]); + else + myParts=particles; + end + nAnal=size(what2Analise,1); + + %% actually analyse + Counts2D=missing(); + CountsX=missing(); CountsY=missing(); + xb=missing(); yb=missing(); + for iAnal=1:nAnal + iColX=strcmpi(declared,what2Analise(iAnal,1)); + if (all(~iColX)), error("...unable to identify coordinate %s!",what2Analise(iAnal,1)); end + iColY=strcmpi(declared,what2Analise(iAnal,2)); + if (all(~iColY)), error("...unable to identify coordinate %s!",what2Analise(iAnal,2)); end + [tmpCounts2D,tmpCountsX,tmpCountsY,tmpXb,tmpYb]=Get2dHistograms(myParts(:,iColX),myParts(:,iColY)); + % store counts + Counts2D=ExpandMat(Counts2D,tmpCounts2D); + CountsX=ExpandMat(CountsX,tmpCountsX'); + CountsY=ExpandMat(CountsY,tmpCountsY'); + xb=ExpandMat(xb,tmpXb'); + yb=ExpandMat(yb,tmpYb'); + end + + %% done + fprintf("...done;\n"); +end diff --git a/optics/BeamStatsShow.m b/optics/BeamStatsShow.m new file mode 100644 index 0000000..2245692 --- /dev/null +++ b/optics/BeamStatsShow.m @@ -0,0 +1,58 @@ +function BeamStatsShow(Counts2D,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle,contours,contourLabels,units,declared) + fprintf("Showing beam statistics...\n"); + + %% set up + if (~exist("lHist","var")), lHist=true; end + if (~exist("myTitle","var")), myTitle=missing(); end + if (~exist("contours","var")), contours=missing(); end + if (~exist("contourLabels","var")), contourLabels=missing(); end + if (~exist("units","var")), units=missing(); end + if (~exist("declared","var")), declared=missing(); end + + if (ismissing(units)), units=["m" "" "m" "" "m" ""]; end + if (ismissing(declared)), declared=["X" "PX" "Y" "PY" "T" "PT"]; end + if (lHist) + % 2D histogram and 1D histograms + ShowMe=Counts2D; + else + % scatter plot and 1D histograms + if (ndims(Counts2D)==3) + ShowMe=reshape(Counts2D,[size(Counts2D,1) size(Counts2D,2)*size(Counts2D,3)]); + else + ShowMe=Counts2D; + end + end + + nAnal=size(what2Analise,1); + + %% actually plot + lSquared=false; + lBinEdges=true; + myConts=missing(); + for iAnal=1:nAnal + iColX=strcmpi(declared,what2Analise(iAnal,1)); + if (all(~iColX)), error("...unable to identify coordinate %s!",what2Analise(iAnal,1)); end + iColY=strcmpi(declared,what2Analise(iAnal,2)); + if (all(~iColY)), error("...unable to identify coordinate %s!",what2Analise(iAnal,2)); end + fullTitle=sprintf("%s vs %s",what2Analise(iAnal,2),what2Analise(iAnal,1)); + if (~ismissing(myTitle)), fullTitle=strcat(myTitle," - ",fullTitle); end + if (~all(ismissing(contours),"all")), myConts=contours(:,:,:,iAnal); end + figure("Name",fullTitle); + xLab=sprintf("%s [%s]",what2Analise(iAnal,1),units(iColX)); + yLab=sprintf("%s [%s]",what2Analise(iAnal,2),units(iColY)); + if (lHist) + % 2D histogram and 1D histograms + [axM,axX,axY]=Plot2DHistograms(ShowMe(:,:,iAnal),CountsX(:,iAnal),CountsY(:,iAnal),xb(:,iAnal),yb(:,iAnal),xLab,yLab,myConts,lHist,lSquared,lBinEdges); + else + % scatter plot and 1D histograms + [axM,axX,axY]=Plot2DHistograms(ShowMe(:,[find(iColX) find(iColY)]),CountsX(:,iAnal),CountsY(:,iAnal),xb(:,iAnal),yb(:,iAnal),xLab,yLab,myConts,lHist,lSquared,lBinEdges); + end + if (~ismissing(contourLabels)) + legend(axM,["Beam" contourLabels],"Location","best"); + end + sgtitle(fullTitle); + end + + %% done + fprintf("...done;\n"); +end From 40fde7f5328c5d81f5d3da06d982723c035f5bd4 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Fri, 7 Jul 2023 16:53:42 +0200 Subject: [PATCH 11/26] Improving numerical format of START condnitions --- MADX-tracking/StartConditionsFormat.m | 16 +++++++++++++--- MADX-tracking/StartConditionsWrite.m | 2 +- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/MADX-tracking/StartConditionsFormat.m b/MADX-tracking/StartConditionsFormat.m index b07f8d7..d20ac00 100644 --- a/MADX-tracking/StartConditionsFormat.m +++ b/MADX-tracking/StartConditionsFormat.m @@ -1,13 +1,23 @@ -function myFmt=StartConditionsFormat(tracker) +function myFmt=StartConditionsFormat(tracker,action) % example format of line of parsed files (argument of ID is read as well): % START, X=-0.01, PX=-0.01, Y=0, PY=0, T=0, PT=0; ! ID=1 % PTC_START, X=-0.01, PX=-0.01, Y=0, PY=0, T=0, PT=0; ! ID=1 + if (~exist("action","var")), action=missing(); end + if (ismissing(action)), action="R"; end switch upper(tracker) case {"MADX","MAD","MAD-X"} - myFmt='START, X=%f, PX=%f, Y=%f, PY=%f, T=%f, PT=%f;%*[^\n]'; + initial="START"; case "PTC" - myFmt='PTC_START, X=%f, PX=%f, Y=%f, PY=%f, T=%f, PT=%f;%*[^\n]'; + initial="PTC_START"; otherwise error("Unknown tracker %s!",tracker); end + switch upper(action) + case {"READ","R"} + myFmt=strcat(initial,", X=%f, PX=%f, Y=%f, PY=%f, T=%f, PT=%f;%*[^\n]"); + case {"WRITE","W"} + myFmt=strcat(initial,", X=% 22.15E, PX=% 22.15E, Y=% 22.15E, PY=% 22.15E, T=% 22.15E, PT=% 22.15E;%*[^\n]"); + otherwise + error("Unknown action %s!",action); + end end diff --git a/MADX-tracking/StartConditionsWrite.m b/MADX-tracking/StartConditionsWrite.m index 54bce51..d1b67f2 100644 --- a/MADX-tracking/StartConditionsWrite.m +++ b/MADX-tracking/StartConditionsWrite.m @@ -1,7 +1,7 @@ function StartConditionsWrite(FileName,particles,tracker) if ( ~exist('tracker','var') ), tracker="MADX"; end fprintf("writing starting conditions of %d particles to file %s (%s format)...\n",length(particles),FileName,tracker); - myFmt=StartConditionsFormat(tracker); + myFmt=StartConditionsFormat(tracker,"W"); fileID = fopen(FileName,'w'); for jj=1:length(particles) fprintf(fileID,"%s\n",sprintf(myFmt,particles(jj,:))); From e676e8d140976e6ac224df533646a62649b93142 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Tue, 1 Aug 2023 16:15:13 +0200 Subject: [PATCH 12/26] a more general way to performing stats on beam samples and showing it --- general/Get2dHistograms.m | 4 +- optics/{BeamContour.m => BeamContourShapes.m} | 4 +- optics/BeamContourStats.m | 21 +++++-- optics/BeamSample.m | 12 ++-- optics/BeamStats.m | 56 +++++++++---------- optics/BeamStatsCalc.m | 35 ++++++++++++ optics/BeamStatsShow.m | 18 ++---- 7 files changed, 94 insertions(+), 56 deletions(-) rename optics/{BeamContour.m => BeamContourShapes.m} (91%) create mode 100644 optics/BeamStatsCalc.m diff --git a/general/Get2dHistograms.m b/general/Get2dHistograms.m index f00f820..f5271ee 100644 --- a/general/Get2dHistograms.m +++ b/general/Get2dHistograms.m @@ -1,6 +1,6 @@ function [nCounts,nCountsX,nCountsY,xb,yb]=Get2dHistograms(currDataX,currDataY,xb,yb,hist1DUsr) - if ( ~exist("xb","var") | ismissing(xb) ), xb=linspace(min(currDataX),max(currDataX),101); end - if ( ~exist("yb","var") | ismissing(yb) ), yb=linspace(min(currDataY),max(currDataY),101); end + if ( ~exist("xb","var") | ismissing(xb) ), xb=linspace(min(currDataX),max(currDataX),101)'; end + if ( ~exist("yb","var") | ismissing(yb) ), yb=linspace(min(currDataY),max(currDataY),101)'; end if ( ~exist("hist1DUsr","var") ), hist1DUsr=true; end nCounts=histcounts2(currDataX,currDataY,xb,yb); if ( hist1DUsr ) diff --git a/optics/BeamContour.m b/optics/BeamContourShapes.m similarity index 91% rename from optics/BeamContour.m rename to optics/BeamContourShapes.m index a06d308..39e695d 100644 --- a/optics/BeamContour.m +++ b/optics/BeamContourShapes.m @@ -1,4 +1,6 @@ -function [Contours]=BeamContour(dType,alfa,beta,emiGeo,iNorm,bb,hh) +function [Contours]=BeamContourShapes(dType,alfa,beta,emiGeo,iNorm,bb,hh) +% BeamContourShapes defines points along defined shapes, for contouring +% beam samples; % % Contours=float(nPoints,2,nPlanes); % NB: iNorm,bb,hh only for bar-of-charge (optional arguments) diff --git a/optics/BeamContourStats.m b/optics/BeamContourStats.m index df3afa9..e98ef7e 100644 --- a/optics/BeamContourStats.m +++ b/optics/BeamContourStats.m @@ -1,9 +1,16 @@ -function [Contours]=BeamContourStats(particles,emiMax) -% particles=float(nParticles,2,nPlanes); +function [Contours]=BeamContourStats(particles,emiMax,declared,what2Analise) +% BeamContourStats defines points along ellypse fitting the beam sample in input +% if emiMax is true, the max single-particle emittance +% is calculated and used instead of the RMS one; +% +% particles=float(nParticles,2*nPlanes); +% + if (~exist("what2Analise","var") | all(ismissing(what2Analise)) ), what2Analise=["X" "PX" ; "Y" "PY" ]; end + if (~exist("declared","var") | all(ismissing(declared)) ), declared=["X" "PX" "Y" "PY" "T" "PT"]; end fprintf("Finding fit ellypse on beam...\n"); %% set up - nPlanes=size(particles,3); + nPlanes=size(what2Analise,2); if (~exist("emiMax","var")), emiMax=false; end if (emiMax) fprintf("...based on largest single-particle geometric emittance!\n"); @@ -14,9 +21,13 @@ %% actually contour Contours=missing(); for iPlane=1:nPlanes - [alpha,beta,emiG]=GetOpticsFromSigmaMatrix(particles(:,:,iPlane)); % ellypse orientation + iColX=find(strcmpi(declared,what2Analise(iPlane,1))); + if (isempty(iColX)), error("...unable to identify coordinate %s!",what2Analise(iPlane,1)); end + iColY=find(strcmpi(declared,what2Analise(iPlane,2))); + if (isempty(iColY)), error("...unable to identify coordinate %s!",what2Analise(iPlane,2)); end + [alpha,beta,emiG]=GetOpticsFromSigmaMatrix(particles(:,[iColX iColY])); % ellypse orientation if (emiMax) - emiG=max(GetSinglePartEmittance(particles(:,:,iPlane),alpha,beta)); % max + emiG=max(GetSinglePartEmittance(particles(:,[iColX iColY]),alpha,beta)); % max end Contours=ExpandMat(Contours,GenPointsAlongEllypse(alpha,beta,emiG)); end diff --git a/optics/BeamSample.m b/optics/BeamSample.m index 566adc8..b072224 100644 --- a/optics/BeamSample.m +++ b/optics/BeamSample.m @@ -1,6 +1,6 @@ function [beamParts]=BeamSample(dType,nParticles,alfa,beta,emiGeo,iNorm,bb,hh) % -% beamParts=float(nParticles,2,nPlanes); +% beamParts=float(nParticles,2*nPlanes); % particle coordinates expressed in physical units, i.e. [m] and [rad]! % NB: iNorm,bb,hh only for bar-of-charge (optional arguments) % @@ -34,19 +34,19 @@ end %% actually sample - beamParts=missing(); iBar=0; + beamParts=NaN(nParticles,2*nPlanes); iBar=0; for iPlane=1:nPlanes switch upper(dType(iPlane)) case {"B","BAR","BAR-OF-CHARGE"} iBar=iBar+1; - beamParts=ExpandMat(beamParts,BeamSampleRect(nParticles,bb(iBar),hh(iBar))); + beamParts(:,2*iPlane-1:2*iPlane)=BeamSampleRect(nParticles,bb(iBar),hh(iBar)); if (iNorm(iBar)) - beamParts(:,:,iPlane)=Norm2Phys(beamParts(:,:,iPlane),beta(iPlane),alfa(iPlane),emiGeo(iPlane)); + beamParts(:,2*iPlane-1:2*iPlane)=Norm2Phys(beamParts(:,2*iPlane-1:2*iPlane),beta(iPlane),alfa(iPlane),emiGeo(iPlane)); end case {"D","DELTA","P","POINT"} - beamParts=ExpandMat(beamParts,zeros(nParticles,2)); + beamParts(:,2*iPlane-1:2*iPlane)=zeros(nParticles,2); case {"G","GAUSS","GAUSSIAN"} - beamParts=ExpandMat(beamParts,BeamSampleGauss(nParticles,alfa(iPlane),beta(iPlane),emiGeo(iPlane))); + beamParts(:,2*iPlane-1:2*iPlane)=BeamSampleGauss(nParticles,alfa(iPlane),beta(iPlane),emiGeo(iPlane)); otherwise error("Unknown distribution type %s!",dType(iPlane)); end diff --git a/optics/BeamStats.m b/optics/BeamStats.m index 7222178..ce9cc3d 100644 --- a/optics/BeamStats.m +++ b/optics/BeamStats.m @@ -1,34 +1,32 @@ -function [Counts2D,CountsX,CountsY,xb,yb]=BeamStats(particles,what2Analise,declared) -% particles=float(nParticles,2,nPlanes); - fprintf("Beam statistics...\n"); +function [Counts2D,CountsX,CountsY,xb,yb]=BeamStats(particles,what2Analise,declared,xbQuery,ybQuery,contours,contourLabels,myTitle) +% +% particles=float(nParticles,2*nPlanes); +% particle coordinates expressed in physical units, i.e. [m] and [rad]! +% + %% interface vars + if (~exist("what2Analise","var") | all(ismissing(what2Analise)) ), what2Analise=["X" "PX" ; "Y" "PY" ]; end + if (~exist("declared","var") | all(ismissing(declared)) ), declared=["X" "PX" "Y" "PY" "T" "PT"]; end + if (~exist("xbQuery","var") | all(ismissing(xbQuery)) ), xbQuery=NaN(1,round(size(particles,2))/2); end % histogram bins set based on data + if (~exist("ybQuery","var") | all(ismissing(ybQuery)) ), ybQuery=NaN(1,round(size(particles,2))/2); end % histogram bins set based on data + if (~exist("contours","var") | ismissing(contours)), contours=missing(); end + if (~exist("contourLabels","var") | ismissing(contourLabels)), contourLabels=missing(); end + if (~exist("myTitle","var") | ismissing(myTitle)), myTitle=missing(); end - %% set up - if (~exist("declared","var")), declared=["X" "PX" "Y" "PY" "T" "PT"]; end - if (ndims(particles)==3) - myParts=reshape(particles,[size(particles,1) size(particles,2)*size(particles,3)]); - else - myParts=particles; - end - nAnal=size(what2Analise,1); + %% contours + contours=ExpandMat(contours,BeamContourStats(particles,false,declared)); % stat ellypses, RMS + contours=ExpandMat(contours,BeamContourStats(particles,true,declared)); % stat ellypses, max single-particle emittance + contours=permute(contours,[1 2 4 3]); + contourLabels=ExpandMat(contourLabels,["RMS" "Max(\epsilon)"],true); % expand array of labels without increasing number of dimensions + + %% actual statistics + [Counts2D,CountsX,CountsY,xb,yb]=BeamStatsCalc(particles,what2Analise,declared,xbQuery,ybQuery); - %% actually analyse - Counts2D=missing(); - CountsX=missing(); CountsY=missing(); - xb=missing(); yb=missing(); - for iAnal=1:nAnal - iColX=strcmpi(declared,what2Analise(iAnal,1)); - if (all(~iColX)), error("...unable to identify coordinate %s!",what2Analise(iAnal,1)); end - iColY=strcmpi(declared,what2Analise(iAnal,2)); - if (all(~iColY)), error("...unable to identify coordinate %s!",what2Analise(iAnal,2)); end - [tmpCounts2D,tmpCountsX,tmpCountsY,tmpXb,tmpYb]=Get2dHistograms(myParts(:,iColX),myParts(:,iColY)); - % store counts - Counts2D=ExpandMat(Counts2D,tmpCounts2D); - CountsX=ExpandMat(CountsX,tmpCountsX'); - CountsY=ExpandMat(CountsY,tmpCountsY'); - xb=ExpandMat(xb,tmpXb'); - yb=ExpandMat(yb,tmpYb'); + %% plot + if (~ismissing(myTitle)) + % 2D histogram + lHist=true; BeamStatsShow(Counts2D,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle,contours,contourLabels,missing(),declared); + % scatter plot + lHist=false; BeamStatsShow(particles,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle,contours,contourLabels,missing(),declared); end - %% done - fprintf("...done;\n"); end diff --git a/optics/BeamStatsCalc.m b/optics/BeamStatsCalc.m new file mode 100644 index 0000000..356b9e0 --- /dev/null +++ b/optics/BeamStatsCalc.m @@ -0,0 +1,35 @@ +function [Counts2D,CountsX,CountsY,xb,yb]=BeamStatsCalc(particles,what2Analise,declared,xbQuery,ybQuery) +% particles=float(nParticles,2*nPlanes); + + fprintf("Beam statistics...\n"); + + %% set up + if (~exist("declared","var") | all(ismissing(declared)) ), declared=["X" "PX" "Y" "PY" "T" "PT"]; end + nAnal=size(what2Analise,1); + % by default, bins of histograms are set based on data + if (~exist("xbQuery","var") | all(ismissing(xbQuery)) ), xbQuery=NaN; end + if (~exist("ybQuery","var") | all(ismissing(ybQuery)) ), ybQuery=NaN; end + + %% actually analyse + Counts2D=missing(); + CountsX=missing(); CountsY=missing(); + xb=missing(); yb=missing(); + for iAnal=1:nAnal + iColX=strcmpi(declared,what2Analise(iAnal,1)); + if (all(~iColX)), error("...unable to identify coordinate %s!",what2Analise(iAnal,1)); end + iColY=strcmpi(declared,what2Analise(iAnal,2)); + if (all(~iColY)), error("...unable to identify coordinate %s!",what2Analise(iAnal,2)); end + if (size(xbQuery,2)==1), xbb=xbQuery; else xbb=xbQuery(:,iAnal); end + if (size(ybQuery,2)==1), ybb=ybQuery; else ybb=ybQuery(:,iAnal); end + [tmpCounts2D,tmpCountsX,tmpCountsY,tmpXb,tmpYb]=Get2dHistograms(particles(:,iColX),particles(:,iColY),xbb,ybb); + % store counts + Counts2D=ExpandMat(Counts2D,tmpCounts2D); + CountsX=ExpandMat(CountsX,tmpCountsX'); + CountsY=ExpandMat(CountsY,tmpCountsY'); + xb=ExpandMat(xb,tmpXb); + yb=ExpandMat(yb,tmpYb); + end + + %% done + fprintf("...done;\n"); +end diff --git a/optics/BeamStatsShow.m b/optics/BeamStatsShow.m index 2245692..35179ed 100644 --- a/optics/BeamStatsShow.m +++ b/optics/BeamStatsShow.m @@ -1,4 +1,7 @@ function BeamStatsShow(Counts2D,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle,contours,contourLabels,units,declared) +% +% Counts2D is either a 2D histogram or a float(nParticles,2*nPlanes) array +% fprintf("Showing beam statistics...\n"); %% set up @@ -11,17 +14,6 @@ function BeamStatsShow(Counts2D,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle if (ismissing(units)), units=["m" "" "m" "" "m" ""]; end if (ismissing(declared)), declared=["X" "PX" "Y" "PY" "T" "PT"]; end - if (lHist) - % 2D histogram and 1D histograms - ShowMe=Counts2D; - else - % scatter plot and 1D histograms - if (ndims(Counts2D)==3) - ShowMe=reshape(Counts2D,[size(Counts2D,1) size(Counts2D,2)*size(Counts2D,3)]); - else - ShowMe=Counts2D; - end - end nAnal=size(what2Analise,1); @@ -42,10 +34,10 @@ function BeamStatsShow(Counts2D,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle yLab=sprintf("%s [%s]",what2Analise(iAnal,2),units(iColY)); if (lHist) % 2D histogram and 1D histograms - [axM,axX,axY]=Plot2DHistograms(ShowMe(:,:,iAnal),CountsX(:,iAnal),CountsY(:,iAnal),xb(:,iAnal),yb(:,iAnal),xLab,yLab,myConts,lHist,lSquared,lBinEdges); + [axM,axX,axY]=Plot2DHistograms(Counts2D(:,:,iAnal),CountsX(:,iAnal),CountsY(:,iAnal),xb(:,iAnal),yb(:,iAnal),xLab,yLab,myConts,lHist,lSquared,lBinEdges); else % scatter plot and 1D histograms - [axM,axX,axY]=Plot2DHistograms(ShowMe(:,[find(iColX) find(iColY)]),CountsX(:,iAnal),CountsY(:,iAnal),xb(:,iAnal),yb(:,iAnal),xLab,yLab,myConts,lHist,lSquared,lBinEdges); + [axM,axX,axY]=Plot2DHistograms(Counts2D(:,[find(iColX) find(iColY)]),CountsX(:,iAnal),CountsY(:,iAnal),xb(:,iAnal),yb(:,iAnal),xLab,yLab,myConts,lHist,lSquared,lBinEdges); end if (~ismissing(contourLabels)) legend(axM,["Beam" contourLabels],"Location","best"); From 63dfca2e2546e9dd83a634229fad367f178568bc Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Tue, 1 Aug 2023 16:15:42 +0200 Subject: [PATCH 13/26] filtering on observation point when parsing dump files by MADX/PTC --- MADX-tracking/ReadSurvivors.m | 79 +++++++++++++++++++---------------- 1 file changed, 44 insertions(+), 35 deletions(-) diff --git a/MADX-tracking/ReadSurvivors.m b/MADX-tracking/ReadSurvivors.m index d1aabe1..bb9eea8 100644 --- a/MADX-tracking/ReadSurvivors.m +++ b/MADX-tracking/ReadSurvivors.m @@ -1,4 +1,4 @@ -function [ particles, obsPoints, nParts ] = ReadSurvivors(path2Files) +function [ particles, obsPoints, nParts ] = ReadSurvivors(path2Files,obsPointNames) % READSURVIVORS Read survivors from MAD-X/PTC tracking. % % example format of line of parsed files (argument of ID is read as well): @@ -14,10 +14,12 @@ % the third dimension identifies the obervation point; % obsPoints: list of names of the obervation points; + %% set up + if (~exist("obsPointNames","var")), obsPointNames=missing(); end [ colNames, colUnits, colFacts, mapping, readFormat ] = ... GetColumnsAndMappingTFS("track"); - % acquire files + %% acquire files particles=missing(); obsPoints=missing(); nParts=missing(); for kk=1:length(path2Files) files=dir(path2Files(kk)); @@ -39,23 +41,27 @@ elseif (startsWith(tLine,"#")) % new section myFields=split(tLine); - nParts=ExpandMat(nParts,double(string(myFields{4}))); - clear tmpParticles; - tmpParticles=cell2mat(textscan(fileID,readFormat)); - fprintf("...acquired %d particles at %s;\n",nParts(end),myFields{end}); - % add closed orbit - if (strcmpi(myFields{end},"START")) - for jj=1:length(CO6Dnames) - if (CO6D(jj)~=0.0) - iCol=mapping(strcmpi(colNames,extractBetween(CO6Dnames(ii),1,strlength(CO6Dnames(ii))-1))); - tmpParticles(:,iCol)=tmpParticles(:,iCol)+CO6D(jj); - fprintf("...taking into account 6D closed orbit on %s at START;\n",colNames(iCol)); + if ( all(ismissing(obsPointNames)) | ~isempty(find(contains(obsPointNames,myFields{end}))) ) + % acquire beam population at the specified obs point + % or at every obs point + nParts=ExpandMat(nParts,double(string(myFields{4}))); + clear tmpParticles; + tmpParticles=cell2mat(textscan(fileID,readFormat)); + fprintf("...acquired %d particles at %s;\n",nParts(end),myFields{end}); + % add closed orbit + if (strcmpi(myFields{end},"START")) + for jj=1:length(CO6Dnames) + if (CO6D(jj)~=0.0) + iCol=mapping(strcmpi(colNames,extractBetween(CO6Dnames(ii),1,strlength(CO6Dnames(ii))-1))); + tmpParticles(:,iCol)=tmpParticles(:,iCol)+CO6D(jj); + fprintf("...taking into account 6D closed orbit on %s at START;\n",colNames(iCol)); + end end end + % store in memory + particles=ExpandMat(particles,tmpParticles); + obsPoints=ExpandMat(obsPoints,string(myFields{end})); end - % store in memory - particles=ExpandMat(particles,tmpParticles); - obsPoints=ExpandMat(obsPoints,string(myFields{end})); end end fclose(fileID); @@ -63,26 +69,29 @@ end fprintf("...done with parsing files;\n"); - % compact storage - uObsPoints=unique(obsPoints); - if (length(uObsPoints)~=length(obsPoints)) - fprintf("...compacting from %d data sets to %d (unique obs points);\n",length(obsPoints),length(uObsPoints)); - readParticles=particles; readObsPoints=obsPoints; readNparts=nParts; - clear particles obsPoints nParts; - particles=missing(); obsPoints=missing(); nParts=missing(); - for ii=1:length(uObsPoints) - iObs=find(strcmp(readObsPoints,uObsPoints(ii))); - clear tmpParticles; - nParts=ExpandMat(nParts,sum(readNparts(iObs))); - tmpParticles=NaN(nParts(end),length(colNames)); - jRead=0; - for jj=1:length(iObs) - tmpParticles(jRead+1:jRead+readNparts(iObs(jj)),:)=readParticles(1:readNparts(iObs(jj)),:,iObs(jj)); - jRead=jRead+readNparts(iObs(jj)); + %% compact storage + % only in case all observation points are parsed + if ( all(ismissing(obsPointNames)) ) + uObsPoints=unique(obsPoints); + if (length(uObsPoints)~=length(obsPoints)) + fprintf("...compacting from %d data sets to %d (unique obs points);\n",length(obsPoints),length(uObsPoints)); + readParticles=particles; readObsPoints=obsPoints; readNparts=nParts; + clear particles obsPoints nParts; + particles=missing(); obsPoints=missing(); nParts=missing(); + for ii=1:length(uObsPoints) + iObs=find(strcmp(readObsPoints,uObsPoints(ii))); + clear tmpParticles; + nParts=ExpandMat(nParts,sum(readNparts(iObs))); + tmpParticles=NaN(nParts(end),length(colNames)); + jRead=0; + for jj=1:length(iObs) + tmpParticles(jRead+1:jRead+readNparts(iObs(jj)),:)=readParticles(1:readNparts(iObs(jj)),:,iObs(jj)); + jRead=jRead+readNparts(iObs(jj)); + end + particles=ExpandMat(particles,tmpParticles); + obsPoints=ExpandMat(obsPoints,uObsPoints(ii)); end - particles=ExpandMat(particles,tmpParticles); - obsPoints=ExpandMat(obsPoints,uObsPoints(ii)); + fprintf("...done with compacting data storage;\n"); end - fprintf("...done with compacting data storage;\n"); end end \ No newline at end of file From 59cc0e25c33cc19c5a060e3a50160458799673a7 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 2 Aug 2023 12:43:35 +0200 Subject: [PATCH 14/26] when computing beam statistics, check that the ellypses (RMS emittance / max single part emittance) are computed only in case of hor/ver phase spaces --- optics/BeamStats.m | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) diff --git a/optics/BeamStats.m b/optics/BeamStats.m index ce9cc3d..946a959 100644 --- a/optics/BeamStats.m +++ b/optics/BeamStats.m @@ -1,4 +1,4 @@ -function [Counts2D,CountsX,CountsY,xb,yb]=BeamStats(particles,what2Analise,declared,xbQuery,ybQuery,contours,contourLabels,myTitle) +function [Counts2D,CountsX,CountsY,xb,yb,contours,contourLabels]=BeamStats(particles,what2Analise,declared,xbQuery,ybQuery,contours,contourLabels,myTitle) % % particles=float(nParticles,2*nPlanes); % particle coordinates expressed in physical units, i.e. [m] and [rad]! @@ -13,10 +13,19 @@ if (~exist("myTitle","var") | ismissing(myTitle)), myTitle=missing(); end %% contours - contours=ExpandMat(contours,BeamContourStats(particles,false,declared)); % stat ellypses, RMS - contours=ExpandMat(contours,BeamContourStats(particles,true,declared)); % stat ellypses, max single-particle emittance - contours=permute(contours,[1 2 4 3]); - contourLabels=ExpandMat(contourLabels,["RMS" "Max(\epsilon)"],true); % expand array of labels without increasing number of dimensions + % stat ellypses, RMS + tmpContour=BeamContourStats(particles,false,declared); + if (any(~ismissing(tmpContour))) + contours=ExpandMat(contours,tmpContour); + contourLabels=ExpandMat(contourLabels,"RMS",true); % expand array of labels without increasing number of dimensions + end + % stat ellypses, max single-part emittance + tmpContour=BeamContourStats(particles,true,declared); + if (any(~ismissing(tmpContour))) + contours=ExpandMat(contours,tmpContour); + contourLabels=ExpandMat(contourLabels,"Max(\epsilon)",true); % expand array of labels without increasing number of dimensions + end + if (any(~ismissing(contours),"all")), contours=permute(contours,[1 2 4 3]); end %% actual statistics [Counts2D,CountsX,CountsY,xb,yb]=BeamStatsCalc(particles,what2Analise,declared,xbQuery,ybQuery); From 80155a744b82c073afbcd8bdf822209fa0e9b7cf Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 2 Aug 2023 12:44:39 +0200 Subject: [PATCH 15/26] forgotten file for previous commit --- optics/BeamContourStats.m | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/optics/BeamContourStats.m b/optics/BeamContourStats.m index e98ef7e..c47883d 100644 --- a/optics/BeamContourStats.m +++ b/optics/BeamContourStats.m @@ -25,11 +25,15 @@ if (isempty(iColX)), error("...unable to identify coordinate %s!",what2Analise(iPlane,1)); end iColY=find(strcmpi(declared,what2Analise(iPlane,2))); if (isempty(iColY)), error("...unable to identify coordinate %s!",what2Analise(iPlane,2)); end - [alpha,beta,emiG]=GetOpticsFromSigmaMatrix(particles(:,[iColX iColY])); % ellypse orientation - if (emiMax) - emiG=max(GetSinglePartEmittance(particles(:,[iColX iColY]),alpha,beta)); % max + if ( any(contains(["X" "Y"],what2Analise(iPlane,1),"IgnoreCase",true)) && any(contains(what2Analise(iPlane,2),"P","IgnoreCase",true)) && ... + (all(contains(what2Analise(iPlane,:),"X","IgnoreCase",true)) || all(contains(what2Analise(iPlane,:),"Y","IgnoreCase",true))) ) + % be sure that we are dealing with hor/ver phase space + [alpha,beta,emiG]=GetOpticsFromSigmaMatrix(particles(:,[iColX iColY])); % ellypse orientation and RMS emittance + if (emiMax), emiG=max(GetSinglePartEmittance(particles(:,[iColX iColY]),alpha,beta)); end % max single-part emittance + Contours=ExpandMat(Contours,GenPointsAlongEllypse(alpha,beta,emiG)); + else + warning("Cannot find an optics ellypse on plane (%s,%s);",what2Analise(iPlane,1),what2Analise(iPlane,2)); end - Contours=ExpandMat(Contours,GenPointsAlongEllypse(alpha,beta,emiG)); end %% done From 3047b8c5786517e7ec5b19bd7d34a82b25c635e0 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 2 Aug 2023 12:45:15 +0200 Subject: [PATCH 16/26] using tiledlayout when plotting 3D spectra --- measurement_analysis/ShowSpectra.m | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/measurement_analysis/ShowSpectra.m b/measurement_analysis/ShowSpectra.m index cc8e850..4f60b92 100644 --- a/measurement_analysis/ShowSpectra.m +++ b/measurement_analysis/ShowSpectra.m @@ -33,18 +33,24 @@ function ShowSpectra(dataSets,tmpTitleFig,addIndex,addLabel,myLabels,myFigSave) nPlanes=length(planes); [nRows,nCols,lDispHor]=GetNrowsNcols(nPlanes*nDataSets,nPlanes); - iPlot=0; + if (lDispHor) + % show planes side by side + tiledlayout(nRows,nCols,'TileSpacing','Compact','Padding','Compact'); % minimise whitespace around plots + else + % show planes on consecutive rows + iPlot=0; + end for iDataSet=1:nDataSets for iPlane=1:nPlanes - iPlot=iPlot+1; if (lDispHor) % show planes side by side - jPlot=iPlot; + nexttile; else % show planes on consecutive rows + iPlot=iPlot+1; jPlot=(iPlane-1)*nCols+(nCols*nPlanes)*(ceil(iPlot/(nCols*nPlanes))-1)+mod(iDataSet-1,nCols)+1; + subplot(nRows,nCols,jPlot); end - subplot(nRows,nCols,jPlot); if ( ~exist('addIndex','var') & ~exist('addLabel','var') ) PlotSpectra(dataSets(:,:,iPlane,iDataSet),BaW); else From ea88bb913848478d6b7282a20a14dfaf5aa15050 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 2 Aug 2023 12:46:00 +0200 Subject: [PATCH 17/26] 3D plot of 2D phase space --- general/Plot2DHistograms.m | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/general/Plot2DHistograms.m b/general/Plot2DHistograms.m index 588c3fd..6bd51f1 100644 --- a/general/Plot2DHistograms.m +++ b/general/Plot2DHistograms.m @@ -21,7 +21,7 @@ axM=subplot('Position', [0.10, 0.10, 0.6, 0.6]); if ( lHist ) hh=histogram2('XBinEdges',xShow,'YBinEdges',yShow,... - 'BinCounts',showMe,'DisplayStyle','tile','ShowEmptyBins','on'); + 'BinCounts',showMe,'FaceColor','flat','ShowEmptyBins','on'); view(2); % starts appearing as a 2D plot else plot(showMe(:,1),showMe(:,2),"k."); @@ -31,7 +31,11 @@ for ii=1:size(contours,3) hold on; contourFmt=".-"; if (ii<=length(contourFmts)), contourFmt=contourFmts(ii); end - plot(contours(:,1,ii),contours(:,2,ii),strcat(contourFmt,"-")); + if (lHist) + plot3(contours(:,1,ii),contours(:,2,ii),ones(size(contours(:,1,ii)))*max(hh.Values,[],"all")/2,strcat(contourFmt,"-")); + else + plot(contours(:,1,ii),contours(:,2,ii),strcat(contourFmt,"-")); + end end end end From 7e39616662350513e0942815ab71ac870896d4be Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 2 Aug 2023 14:08:34 +0200 Subject: [PATCH 18/26] adding a function for comparing spectra belonging to different scans one-by-one --- measurement_analysis/ShowSpectraCompare.m | 74 +++++++++++++++++++++++ 1 file changed, 74 insertions(+) create mode 100644 measurement_analysis/ShowSpectraCompare.m diff --git a/measurement_analysis/ShowSpectraCompare.m b/measurement_analysis/ShowSpectraCompare.m new file mode 100644 index 0000000..f86c772 --- /dev/null +++ b/measurement_analysis/ShowSpectraCompare.m @@ -0,0 +1,74 @@ +function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg) +% profiles=float(nBinCentres,1+nProfiles,nSeries) + + %% interface vars + if ( ~exist("manipulate","var") | all(ismissing(manipulate)) ), manipulate=missing(); end + if ( ~exist("yLab","var") | ismissing(yLab) ), yLab=missing(); end + if ( ~exist("xLab","var") | ismissing(xLab) ), xLab=missing(); end + if ( ~exist("myLeg","var") | ismissing(myLeg) ), myLeg=missing(); end + if (~ismissing(manipulate)) + for iMan=1:length(manipulate) + switch upper(manipulate(iMan)) + case {"CENTRE","CENTER"} + if ( ismissing(xLab) ), xLab="[mm]"; end + case {"AREA","ONE"} + if ( ismissing(yLab) ), yLab="pdf (Area=1) []"; end + case {"PEAK","MAX"} + if ( ismissing(yLab) ), yLab="pdf (norm to peaks) []"; end + otherwise + error("unknown manipulation: %s!",manipulate(iMan)); + end + end + end + if ( ismissing(yLab) ), yLab="counts []"; end + if ( ismissing(xLab) ), xLab="[mm]"; end + + %% set up + nProfiles=size(profiles,2)-1; + nDataSets=size(profiles,3); + if (~ismissing(myLeg)) + [nRows,nCols]=GetNrowsNcols(nProfiles+1); + else + [nRows,nCols]=GetNrowsNcols(nProfiles); + end + ff=figure(); + ff.Position(1:2)=[0 0]; % figure at lower-left corner of screen + ff.Position(3)=ff.Position(3)*nCols/2; % larger figure + ff.Position(4)=ff.Position(4)*nRows/2; % larger figure + cm=colormap(jet(nDataSets)); + tiledlayout(nRows,nCols,'TileSpacing','Compact','Padding','Compact'); % minimise whitespace around plots + + %% actually plot + % - profiles + for iProfile=1:nProfiles + nexttile; + for iDataSet=1:nDataSets + if (iDataSet>1), hold on; end + showX=profiles(:,1,iDataSet); + showY=profiles(:,1+iProfile,iDataSet); + if (~ismissing(manipulate)) + for iMan=1:length(manipulate) + switch upper(manipulate(iMan)) + case {"CENTRE","CENTER"} + showX=showX-mean(showX); + case {"AREA","ONE"} + showY=showY/sum(showY); + case {"PEAK","MAX"} + showY=showY/max(showY); + end + end + end + bar(showX,showY,1,"FaceColor","none","EdgeColor",cm(iDataSet,:)); + end + xlabel(xLab); ylabel(yLab); grid on; title(sprintf("profile %d",iProfile)); + end + % - legend plot + if (~ismissing(myLeg)) + nexttile; + for iDataSet=1:nDataSets + if (iDataSet>1), hold on; end + plot(NaN(),NaN(),".-","Color",cm(iDataSet,:)); hold on; + end + legend(myLeg,"Location","best"); + end +end From 26f317f1c019e6a8f14c15956a25bf1b180a71f5 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 2 Aug 2023 16:52:13 +0200 Subject: [PATCH 19/26] different way of indexing dimensions in BeamSample.m --- optics/BeamSample.m | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/optics/BeamSample.m b/optics/BeamSample.m index b072224..70577b0 100644 --- a/optics/BeamSample.m +++ b/optics/BeamSample.m @@ -39,14 +39,14 @@ switch upper(dType(iPlane)) case {"B","BAR","BAR-OF-CHARGE"} iBar=iBar+1; - beamParts(:,2*iPlane-1:2*iPlane)=BeamSampleRect(nParticles,bb(iBar),hh(iBar)); + beamParts(:,(1:2)+2*(iPlane-1))=BeamSampleRect(nParticles,bb(iBar),hh(iBar)); if (iNorm(iBar)) - beamParts(:,2*iPlane-1:2*iPlane)=Norm2Phys(beamParts(:,2*iPlane-1:2*iPlane),beta(iPlane),alfa(iPlane),emiGeo(iPlane)); + beamParts(:,(1:2)+2*(iPlane-1))=Norm2Phys(beamParts(:,(1:2)+2*(iPlane-1)),beta(iPlane),alfa(iPlane),emiGeo(iPlane)); end case {"D","DELTA","P","POINT"} - beamParts(:,2*iPlane-1:2*iPlane)=zeros(nParticles,2); + beamParts(:,(1:2)+2*(iPlane-1))=zeros(nParticles,2); case {"G","GAUSS","GAUSSIAN"} - beamParts(:,2*iPlane-1:2*iPlane)=BeamSampleGauss(nParticles,alfa(iPlane),beta(iPlane),emiGeo(iPlane)); + beamParts(:,(1:2)+2*(iPlane-1))=BeamSampleGauss(nParticles,alfa(iPlane),beta(iPlane),emiGeo(iPlane)); otherwise error("Unknown distribution type %s!",dType(iPlane)); end From 6051b9f5a9ca2f1e039b11293fdb50ff82262199 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 2 Aug 2023 16:52:58 +0200 Subject: [PATCH 20/26] minor changes in ShowSpectraCompare.m: * autoranging x-axis based on the leftmost and rightmost non-zero bin; * decreasing thickness of histogram in each plot; * dedicated plot with legend only created only if we have more than 6 plots; * plot titles can be passed through the function interface; --- measurement_analysis/ShowSpectraCompare.m | 21 ++++++++++++++++----- 1 file changed, 16 insertions(+), 5 deletions(-) diff --git a/measurement_analysis/ShowSpectraCompare.m b/measurement_analysis/ShowSpectraCompare.m index f86c772..11d70d9 100644 --- a/measurement_analysis/ShowSpectraCompare.m +++ b/measurement_analysis/ShowSpectraCompare.m @@ -1,4 +1,4 @@ -function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg) +function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg,myProfLabels) % profiles=float(nBinCentres,1+nProfiles,nSeries) %% interface vars @@ -6,6 +6,7 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg) if ( ~exist("yLab","var") | ismissing(yLab) ), yLab=missing(); end if ( ~exist("xLab","var") | ismissing(xLab) ), xLab=missing(); end if ( ~exist("myLeg","var") | ismissing(myLeg) ), myLeg=missing(); end + if ( ~exist("myProfLabels","var") | all(ismissing(myProfLabels)) ), myProfLabels=missing(); end if (~ismissing(manipulate)) for iMan=1:length(manipulate) switch upper(manipulate(iMan)) @@ -26,7 +27,7 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg) %% set up nProfiles=size(profiles,2)-1; nDataSets=size(profiles,3); - if (~ismissing(myLeg)) + if (nProfiles>6) [nRows,nCols]=GetNrowsNcols(nProfiles+1); else [nRows,nCols]=GetNrowsNcols(nProfiles); @@ -42,6 +43,7 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg) % - profiles for iProfile=1:nProfiles nexttile; + myRange=NaN(nDataSets,2); for iDataSet=1:nDataSets if (iDataSet>1), hold on; end showX=profiles(:,1,iDataSet); @@ -58,12 +60,21 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg) end end end - bar(showX,showY,1,"FaceColor","none","EdgeColor",cm(iDataSet,:)); + bar(showX,showY,1,"FaceColor","none","EdgeColor",cm(iDataSet,:),"LineWidth",1/iDataSet); + myRange(iDataSet,1)=showX(find(showY>0,1,"first")); % min xVal with a non-zero yVal + myRange(iDataSet,2)=showX(find(showY>0,1,"last")); % max xVal with a non-zero yVal end - xlabel(xLab); ylabel(yLab); grid on; title(sprintf("profile %d",iProfile)); + xlabel(xLab); ylabel(yLab); grid on; + if (any(~ismissing(myProfLabels))) + title(myProfLabels(iProfile)); + else + title(sprintf("profile %d",iProfile)); + end + xlim([min(myRange(:,1)) max(myRange(:,2))]); + if (nProfiles<=6), legend(myLeg,"Location","best"); end end % - legend plot - if (~ismissing(myLeg)) + if (nProfiles>6) nexttile; for iDataSet=1:nDataSets if (iDataSet>1), hold on; end From dda266196b967afd32e3acc9b8e0d2a0e88a2b78 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 3 Aug 2023 16:07:34 +0200 Subject: [PATCH 21/26] adding sinogram-like possibility to plot spectra --- measurement_analysis/PlotSinogram.m | 34 +++++++++++++++++ measurement_analysis/PlotSpectra.m | 57 +++++++++++++++++++---------- measurement_analysis/ShowSpectra.m | 36 +++++++++++++----- 3 files changed, 99 insertions(+), 28 deletions(-) create mode 100644 measurement_analysis/PlotSinogram.m diff --git a/measurement_analysis/PlotSinogram.m b/measurement_analysis/PlotSinogram.m new file mode 100644 index 0000000..3fd097b --- /dev/null +++ b/measurement_analysis/PlotSinogram.m @@ -0,0 +1,34 @@ +function PlotSinogram(xs,ys,data,l3D) +% input variables: +% - xs, ys: bin centres (1D float); +% - data: 2D (color) map (2D float); +% . data vs X: cols; +% . data vs Y: rows; +% - l3D: boolean requesting a 3D plot (top view by default); + if (~exist("l3D","var")), l3D=true; end + [xT,idx]=sort(xs); + [yT,idy]=sort(ys); + if (l3D) + % 3D plot via histogram2: + % - XBinEdges,YBinEdges are bin-edge values! + % - size(BinCounts)=[length(XBinEdges)-1 length(YBinEdges)-1] + xShow=BinCentres2Edges(xT); + yShow=BinCentres2Edges(yT); + showMe=data'; + showMe=showMe(idx,:); + showMe=showMe(:,idy); + histogram2('XBinEdges',xShow,'YBinEdges',yShow,... + 'BinCounts',showMe,'FaceColor','flat','ShowEmptyBins','on'); + view(2); % starts appearing as a 2D plot + else + % 2D plot via imagesc: + % - XData,YData are centre values! + % - size(CData)=[length(YData) length(XData)] + showMe=data; + showMe=showMe(:,idx); + showMe=showMe(idy,:); + imagesc('XData',xT,'YData',yT,'CData',showMe); + colorbar; + end + grid(); +end diff --git a/measurement_analysis/PlotSpectra.m b/measurement_analysis/PlotSpectra.m index c5940f9..f333b0b 100644 --- a/measurement_analysis/PlotSpectra.m +++ b/measurement_analysis/PlotSpectra.m @@ -1,4 +1,4 @@ -function PlotSpectra(dataSets,BaW,addIndex,addLabel) +function PlotSpectra(dataSets,BaW,addIndex,addLabel,iMod) % PlotSpectra plots distributions/histograms in a 3D space with a free parameter; % this plotting function can be useful to explore eg % spectra taken during scans; @@ -15,29 +15,48 @@ function PlotSpectra(dataSets,BaW,addIndex,addLabel) % - addIndex [float(nColumns-1), optional]: list of IDs to be shown; % it can be used to separate distribution by cyProg or cyCode; % - addLabel [string, optional]: name of the y-axis; +% - iMod [integer, optional]: type of plot: +% . 1: colored histograms in a 3D view; +% . 2: sinogram-like view; +% . 3: 3D sinogram-like; % % see also ParseSFMData, ShowSpectra and SumSpectra. % + + %% pre-processing nDataSets=size(dataSets,2)-1; - if ( ~exist('BaW','var') ) - BaW=false; - end - if ( BaW ) - cm=gray(nDataSets); - else - cm=colormap(parula(nDataSets)); - end - if ( ~exist('addIndex','var') ) - addIndex=1:nDataSets; - end - if ( ~exist('addLabel','var') ) - addLabel="ID"; - end - for iSet=1:nDataSets - PlotSpectrum(dataSets(:,1),dataSets(:,1+iSet),addIndex(iSet),cm(iSet,:)); - hold on; + if ( ~exist('BaW','var') ), BaW=false; end + if ( ~exist('addIndex','var') | all(ismissing(addIndex)) ), addIndex=1:nDataSets; end + if ( ~exist('addLabel','var') | all(ismissing(addLabel)) ), addLabel="ID"; end + if ( ~exist('iMod','var') | ismissing(iMod) ), iMod=1; end + + %% actually plot + switch iMod + case 1 + % colored histograms in a 3D view + if ( BaW ) + cm=gray(nDataSets); + else + cm=colormap(parula(nDataSets)); + end + for iSet=1:nDataSets + if (iSet>1), hold on; end + PlotSpectrum(dataSets(:,1),dataSets(:,1+iSet),addIndex(iSet),cm(iSet,:)); + end + ylabel(LabelMe(addLabel)); + case 2 + % sinogram-like view + l3D=false; + PlotSinogram(addIndex,dataSets(:,1),dataSets(:,2:end),l3D); + xlabel(LabelMe(addLabel)); + case 3 + % 3D sinogram-like + l3D=true; + PlotSinogram(addIndex,dataSets(:,1),dataSets(:,2:end),l3D); + xlabel(LabelMe(addLabel)); + otherwise + error("Wrong mode! %d",iMod); end - ylabel(LabelMe(addLabel)); grid on; end diff --git a/measurement_analysis/ShowSpectra.m b/measurement_analysis/ShowSpectra.m index 4f60b92..5e2233b 100644 --- a/measurement_analysis/ShowSpectra.m +++ b/measurement_analysis/ShowSpectra.m @@ -1,4 +1,4 @@ -function ShowSpectra(dataSets,tmpTitleFig,addIndex,addLabel,myLabels,myFigSave) +function ShowSpectra(dataSets,tmpTitleFig,addIndex,addLabel,myLabels,myFigSave,iMod) % ShowSpectra shows distributions recorded by SFM, QBM and GIM; % it shows a 1x2 3D figure, with the distributions on the % horizontal plane on the left and those on the vertical @@ -21,11 +21,23 @@ function ShowSpectra(dataSets,tmpTitleFig,addIndex,addLabel,myLabels,myFigSave) % - addLabel [string, optional]: name of the y-axis; % - myLabels [strings(nSets), optional]: a label for each scan; % - myFigSave (string, optional): file name where to save the plot; +% - iMod [integer, optional]: type of plot: +% . 1: colored histograms in a 3D view; +% . 2: sinogram-like view; +% . 3: 3D sinogram-like; % % see also ParseSFMData, PlotSpectra and SumSpectra. fprintf("plotting data...\n"); + + %% pre-processing + if (~exist("addIndex","var") | all(ismissing(addIndex))), addIndex=missing(); end + if (~exist("addLabel","var") | all(ismissing(addLabel))), addLabel=missing(); end + if (~exist("myLabels","var") | all(ismissing(addLabel))), myLabels=missing(); end if (~exist("myFigSave","var")), myFigSave=missing(); end + if (~exist('iMod','var') | ismissing(iMod)), iMod=1; end + + %% actually plotting ff=figure('Name',LabelMe(tmpTitleFig),'NumberTitle','off'); BaW=false; % always use colored plots nDataSets=size(dataSets,4); @@ -51,18 +63,24 @@ function ShowSpectra(dataSets,tmpTitleFig,addIndex,addLabel,myLabels,myFigSave) jPlot=(iPlane-1)*nCols+(nCols*nPlanes)*(ceil(iPlot/(nCols*nPlanes))-1)+mod(iDataSet-1,nCols)+1; subplot(nRows,nCols,jPlot); end - if ( ~exist('addIndex','var') & ~exist('addLabel','var') ) - PlotSpectra(dataSets(:,:,iPlane,iDataSet),BaW); + if ( ismissing(addIndex) & ismissing(addLabel) ) + PlotSpectra(dataSets(:,:,iPlane,iDataSet),BaW,missing(),missing(),iMod); else - PlotSpectra(dataSets(:,:,iPlane,iDataSet),BaW,addIndex(:,iDataSet),addLabel); + PlotSpectra(dataSets(:,:,iPlane,iDataSet),BaW,addIndex(:,iDataSet),addLabel,iMod); end - if (exist('myLabels','var')) - title(sprintf("%s - %s",myLabels(iDataSet),planes(iPlane))); - else + if (ismissing(myLabels)) title(planes(iPlane)); + else + title(sprintf("%s - %s",myLabels(iDataSet),planes(iPlane))); + end + switch iMod + case 1 + xlabel("position [mm]"); + zlabel("Counts []"); + otherwise + ylabel("position [mm]"); + zlabel("Counts []"); end - xlabel("position [mm]"); - zlabel("Counts []"); end end From 0aea23870b1ca00a194756533cec476222065ee8 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 3 Aug 2023 16:08:25 +0200 Subject: [PATCH 22/26] bug fix in ShowSpectraCompare.m when looking for x-range + better legend + first data set is taken as ref in black solid --- measurement_analysis/ShowSpectraCompare.m | 24 ++++++++++++++--------- 1 file changed, 15 insertions(+), 9 deletions(-) diff --git a/measurement_analysis/ShowSpectraCompare.m b/measurement_analysis/ShowSpectraCompare.m index 11d70d9..10bf522 100644 --- a/measurement_analysis/ShowSpectraCompare.m +++ b/measurement_analysis/ShowSpectraCompare.m @@ -26,8 +26,8 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg,myProfLabels) %% set up nProfiles=size(profiles,2)-1; - nDataSets=size(profiles,3); - if (nProfiles>6) + nDataSets=size(profiles,3)-1; % first profile taken as "reference" + if (nDataSets>3) [nRows,nCols]=GetNrowsNcols(nProfiles+1); else [nRows,nCols]=GetNrowsNcols(nProfiles); @@ -60,9 +60,15 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg,myProfLabels) end end end - bar(showX,showY,1,"FaceColor","none","EdgeColor",cm(iDataSet,:),"LineWidth",1/iDataSet); - myRange(iDataSet,1)=showX(find(showY>0,1,"first")); % min xVal with a non-zero yVal - myRange(iDataSet,2)=showX(find(showY>0,1,"last")); % max xVal with a non-zero yVal + if (iDataSet==1) + bar(showX,showY,1,"FaceColor","black","EdgeColor","black"); + else + bar(showX,showY,1,"FaceColor","none","EdgeColor",cm(iDataSet-1,:),"LineWidth",1/iDataSet); + end + myMin=find(showY>0,1,"first"); % min xVal with a non-zero yVal + if (isempty(myMin)), myRange(iDataSet,1)=min(showX); else, myRange(iDataSet,1)=showX(myMin); end + myMax=find(showY>0,1,"last"); % max xVal with a non-zero yVal + if (isempty(myMax)), myRange(iDataSet,2)=max(showX); else, myRange(iDataSet,2)=showX(myMax); end end xlabel(xLab); ylabel(yLab); grid on; if (any(~ismissing(myProfLabels))) @@ -71,14 +77,14 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg,myProfLabels) title(sprintf("profile %d",iProfile)); end xlim([min(myRange(:,1)) max(myRange(:,2))]); - if (nProfiles<=6), legend(myLeg,"Location","best"); end + if (nDataSets<=3 & ~ismissing(myLeg) & iProfile==1), legend(myLeg,"Location","best"); end end % - legend plot - if (nProfiles>6) + if (nDataSets>3 & ~ismissing(myLeg)) nexttile; + bar(NaN(),NaN(),1,"FaceColor","black"); for iDataSet=1:nDataSets - if (iDataSet>1), hold on; end - plot(NaN(),NaN(),".-","Color",cm(iDataSet,:)); hold on; + hold on; bar(NaN(),NaN(),1,"FaceColor","none","EdgeColor",cm(iDataSet,:)); end legend(myLeg,"Location","best"); end From 90a2cf8bf5c61cee65c25a647caeb2259aafd2bb Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 3 Aug 2023 16:09:45 +0200 Subject: [PATCH 23/26] units can be passed through interface of BeamStats.m function --- optics/BeamStats.m | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/optics/BeamStats.m b/optics/BeamStats.m index 946a959..c288a72 100644 --- a/optics/BeamStats.m +++ b/optics/BeamStats.m @@ -1,4 +1,4 @@ -function [Counts2D,CountsX,CountsY,xb,yb,contours,contourLabels]=BeamStats(particles,what2Analise,declared,xbQuery,ybQuery,contours,contourLabels,myTitle) +function [Counts2D,CountsX,CountsY,xb,yb,contours,contourLabels]=BeamStats(particles,what2Analise,declared,units,xbQuery,ybQuery,contours,contourLabels,myTitle) % % particles=float(nParticles,2*nPlanes); % particle coordinates expressed in physical units, i.e. [m] and [rad]! @@ -6,6 +6,7 @@ %% interface vars if (~exist("what2Analise","var") | all(ismissing(what2Analise)) ), what2Analise=["X" "PX" ; "Y" "PY" ]; end if (~exist("declared","var") | all(ismissing(declared)) ), declared=["X" "PX" "Y" "PY" "T" "PT"]; end + if (~exist("units","var") | all(ismissing(units))), units=["m" "" "m" "" "m" ""]; end if (~exist("xbQuery","var") | all(ismissing(xbQuery)) ), xbQuery=NaN(1,round(size(particles,2))/2); end % histogram bins set based on data if (~exist("ybQuery","var") | all(ismissing(ybQuery)) ), ybQuery=NaN(1,round(size(particles,2))/2); end % histogram bins set based on data if (~exist("contours","var") | ismissing(contours)), contours=missing(); end @@ -33,9 +34,9 @@ %% plot if (~ismissing(myTitle)) % 2D histogram - lHist=true; BeamStatsShow(Counts2D,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle,contours,contourLabels,missing(),declared); + lHist=true; BeamStatsShow(Counts2D,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle,contours,contourLabels,units,declared); % scatter plot - lHist=false; BeamStatsShow(particles,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle,contours,contourLabels,missing(),declared); + lHist=false; BeamStatsShow(particles,CountsX,CountsY,xb,yb,what2Analise,lHist,myTitle,contours,contourLabels,units,declared); end end From d20439974d47e82103241c229a0f0e5854e98ad3 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 24 Aug 2023 10:09:39 +0200 Subject: [PATCH 24/26] considering also CO in stat analisis of a beam sample --- optics/BeamContourShapes.m | 23 +++++++++++++++++++++-- optics/BeamContourStats.m | 6 +++--- optics/BeamSample.m | 18 +++++++++++++++++- optics/GetOpticsFromSigmaMatrix.m | 4 +++- 4 files changed, 44 insertions(+), 7 deletions(-) diff --git a/optics/BeamContourShapes.m b/optics/BeamContourShapes.m index 39e695d..7719f30 100644 --- a/optics/BeamContourShapes.m +++ b/optics/BeamContourShapes.m @@ -1,8 +1,9 @@ -function [Contours]=BeamContourShapes(dType,alfa,beta,emiGeo,iNorm,bb,hh) +function [Contours]=BeamContourShapes(dType,alfa,beta,emiGeo,CO,iNorm,bb,hh) % BeamContourShapes defines points along defined shapes, for contouring % beam samples; % % Contours=float(nPoints,2,nPlanes); +% CO added only if provided % NB: iNorm,bb,hh only for bar-of-charge (optional arguments) % fprintf("Contouring beam...\n"); @@ -19,6 +20,15 @@ if (nPlanes~=length(emiGeo)) error("Incosistent number of types (%d) and geometric emittances (%d) (including NaNs)!",nPlanes,length(emiGeo)); end + if (~exist("CO","var") | all(ismissing(CO))) + CO=zeros(2*nPlanes,1); + else + if (length(CO)~=2*nPlanes) + error("Incosistent number of dimensions (%d) and CO specifications (%d)!",2*nPlanes,length(CO)); + end + % avoid NaNs + CO(ismissing(CO))=0.0; + end % - specific to BAR-of-charge nBars=sum(strcmpi(dType,"BAR")); if (~exist("iNorm","var")), iNorm=false(nBars,1); end @@ -35,10 +45,11 @@ end %% actually contour - Contours=missing(); + Contours=missing(); iBar=0; for iPlane=1:nPlanes switch upper(dType(iPlane)) case {"B","BAR","BAR-OF-CHARGE"} + iBar=iBar+1; Contours=ExpandMat(Contours,GenPointsAlongRectangle(bb(iBar),hh(iBar))); if (iNorm(iBar)) Contours(:,:,iPlane)=Norm2Phys(Contours(:,:,iPlane),beta(iPlane),alfa(iPlane),emiGeo(iPlane)); @@ -52,6 +63,14 @@ end end + %% add closed orbit + if (any(CO~=0.0)) + fprintf("...adding closed orbit...;\n"); + for iPlane=1:nPlanes + Contours(:,:,iPlane)=Contours(:,:,iPlane)+CO((1:2)+2*(iPlane-1)); + end + end + %% done fprintf("...done;\n"); end diff --git a/optics/BeamContourStats.m b/optics/BeamContourStats.m index c47883d..eb52550 100644 --- a/optics/BeamContourStats.m +++ b/optics/BeamContourStats.m @@ -28,9 +28,9 @@ if ( any(contains(["X" "Y"],what2Analise(iPlane,1),"IgnoreCase",true)) && any(contains(what2Analise(iPlane,2),"P","IgnoreCase",true)) && ... (all(contains(what2Analise(iPlane,:),"X","IgnoreCase",true)) || all(contains(what2Analise(iPlane,:),"Y","IgnoreCase",true))) ) % be sure that we are dealing with hor/ver phase space - [alpha,beta,emiG]=GetOpticsFromSigmaMatrix(particles(:,[iColX iColY])); % ellypse orientation and RMS emittance - if (emiMax), emiG=max(GetSinglePartEmittance(particles(:,[iColX iColY]),alpha,beta)); end % max single-part emittance - Contours=ExpandMat(Contours,GenPointsAlongEllypse(alpha,beta,emiG)); + [alpha,beta,emiG,~,bars]=GetOpticsFromSigmaMatrix(particles(:,[iColX iColY])); % ellypse orientation and RMS emittance + if (emiMax), emiG=max(GetSinglePartEmittance(particles(:,[iColX iColY])-bars,alpha,beta)); end % max single-part emittance + Contours=ExpandMat(Contours,GenPointsAlongEllypse(alpha,beta,emiG)+bars); else warning("Cannot find an optics ellypse on plane (%s,%s);",what2Analise(iPlane,1),what2Analise(iPlane,2)); end diff --git a/optics/BeamSample.m b/optics/BeamSample.m index 70577b0..9c6e8dc 100644 --- a/optics/BeamSample.m +++ b/optics/BeamSample.m @@ -1,7 +1,8 @@ -function [beamParts]=BeamSample(dType,nParticles,alfa,beta,emiGeo,iNorm,bb,hh) +function [beamParts]=BeamSample(dType,nParticles,alfa,beta,emiGeo,CO,iNorm,bb,hh) % % beamParts=float(nParticles,2*nPlanes); % particle coordinates expressed in physical units, i.e. [m] and [rad]! +% CO added only if provided % NB: iNorm,bb,hh only for bar-of-charge (optional arguments) % fprintf("Generating beam of %d particles...\n",nParticles); @@ -18,6 +19,15 @@ if (nPlanes~=length(emiGeo)) error("Incosistent number of types (%d) and geometric emittances (%d) (including NaNs)!",nPlanes,length(emiGeo)); end + if (~exist("CO","var") | all(ismissing(CO))) + CO=zeros(2*nPlanes,1); + else + if (length(CO)~=2*nPlanes) + error("Incosistent number of dimensions (%d) and CO specifications (%d)!",2*nPlanes,length(CO)); + end + % avoid NaNs + CO(ismissing(CO))=0.0; + end % - specific to BAR-of-charge nBars=sum(strcmpi(dType,"BAR")); if (~exist("iNorm","var")), iNorm=false(nBars,1); end @@ -52,6 +62,12 @@ end end + %% add closed orbit + if (any(CO~=0.0)) + fprintf("...adding closed orbit...;\n"); + beamParts=beamParts+CO; + end + %% done fprintf("...done;\n"); end diff --git a/optics/GetOpticsFromSigmaMatrix.m b/optics/GetOpticsFromSigmaMatrix.m index 2214999..407a779 100644 --- a/optics/GetOpticsFromSigmaMatrix.m +++ b/optics/GetOpticsFromSigmaMatrix.m @@ -1,4 +1,4 @@ -function [alpha,beta,emiG,sigM]=GetOpticsFromSigmaMatrix(MyPoints) +function [alpha,beta,emiG,sigM,bars]=GetOpticsFromSigmaMatrix(MyPoints) % GetOpticsFromSigmaMatrix(MyPoints) to get the beam optics functions from % the statistical analysis of the beam % @@ -8,8 +8,10 @@ % - alpha,beta,emiG (float): optics functions and emittance of the beam % population [,m,m rad]; % - sigM (float(2,2)): sigma matrix; +% - bars (float(2,1)): mean of the point distributions; % one plane at a time! + bars=mean(MyPoints); sigM=cov(MyPoints); [beta,alpha,emiG]=DecodeOpticsFit([sigM(1,1) sigM(1,2) sigM(2,2)]); end From ee8d488e808df7d763662d23d2bf0d86e8ff943f Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 24 Aug 2023 10:10:38 +0200 Subject: [PATCH 25/26] correct number of header lines when parsing a .tfs table created by PTC --- MADX-optics/ParseTfsTable.m | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/MADX-optics/ParseTfsTable.m b/MADX-optics/ParseTfsTable.m index 9015985..51720be 100644 --- a/MADX-optics/ParseTfsTable.m +++ b/MADX-optics/ParseTfsTable.m @@ -41,7 +41,9 @@ else switch upper(genBy) case "TWISS" - nHeader=48; + nHeader=48; % 52; % 48; + case "PTC_TWISS" + nHeader=90; case "SCAN" nHeader=1; otherwise From 09603b105a42195896b2517b483def9d1acb768e Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 24 Aug 2023 10:11:37 +0200 Subject: [PATCH 26/26] adding some printout on terminal when comparing spectra --- measurement_analysis/ShowSpectraCompare.m | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/measurement_analysis/ShowSpectraCompare.m b/measurement_analysis/ShowSpectraCompare.m index 10bf522..c040a12 100644 --- a/measurement_analysis/ShowSpectraCompare.m +++ b/measurement_analysis/ShowSpectraCompare.m @@ -39,6 +39,8 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg,myProfLabels) cm=colormap(jet(nDataSets)); tiledlayout(nRows,nCols,'TileSpacing','Compact','Padding','Compact'); % minimise whitespace around plots + fprintf("Comparing %d profiles from %d datasets...\n",nProfiles,nDataSets); + %% actually plot % - profiles for iProfile=1:nProfiles @@ -88,4 +90,6 @@ function ShowSpectraCompare(profiles,manipulate,xLab,yLab,myLeg,myProfLabels) end legend(myLeg,"Location","best"); end + + fprintf("...done;\n"); end