From 2e34400087200909b132d285357b8f497fc5b7e1 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 31 Aug 2023 17:10:40 +0200 Subject: [PATCH 1/4] `GetColumnsAndMappingTFS.m` now also handles start contiditions of PTC/MADX tracking --- MADX-optics/GetColumnsAndMappingTFS.m | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/MADX-optics/GetColumnsAndMappingTFS.m b/MADX-optics/GetColumnsAndMappingTFS.m index 90a91a5..22c5fca 100644 --- a/MADX-optics/GetColumnsAndMappingTFS.m +++ b/MADX-optics/GetColumnsAndMappingTFS.m @@ -116,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','RE'} + case {'RM','RMATRIX','RE','TM','TMATRIX'} colNames=[ "NAME" "KEYWORD" "L" "S" ... "RE11" "RE12" "RE21" "RE22" "RE16" "RE26" ... "RE33" "RE34" "RE43" "RE44" "RE36" "RE46" ... @@ -142,6 +142,10 @@ "m" "" "m" "" "m" "" ... "m" "GeV" ]; readFormat = '%f %f %f %f %f %f %f %f %f %f'; + case {'START','STARTING'} + colNames=[ "X" "PX" "Y" "PY" "T" "PT" ]; + colUnits=[ "m" "" "m" "" "m" "" ]; + readFormat = '%f %f %f %f %f %f'; otherwise error('which column mapping for TFS table?'); end @@ -164,7 +168,7 @@ '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... '%f',mySep,'%f',mySep,'%f',mySep,'%f',mySep, ... '%f',mySep,'%f' ); - case {'GEO','GEOMETRY'} + case {'GEO','GEOMETRY','GEOM'} colNames=[ "BRHO" "BP" "ID" ... "KICK" "HKICK" "VKICK" "ANGLE" "K0L" "K1L" "K2L" ... "APERTYPE" "APER_1" "APER_2" "APER_3" "APER_4" "APOFF_1" "APOFF_2" ]; @@ -174,7 +178,7 @@ 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'} + case {'RM','RMATRIX','RE','TM','TMATRIX'} colNames=[ "BRHO" "BP" "ID" ... "RE11" "RE12" "RE21" "RE22" "RE16" "RE26" ... "RE33" "RE34" "RE43" "RE44" "RE36" "RE46" ... @@ -187,7 +191,7 @@ '%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'} + case {'LOSS','LOSSES','LOST','TRACK','TRACKS','TRACKING','SURV','SURVS','SURVIVOR','SURVIVORS','START','STARTING'} error('... %s .tfs type generated by %s NOT available yet!',whichTFS,genBy); otherwise error('which column mapping for TFS table?'); From 16142825479aff77d8ec001fd03516f4bc03a878 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Fri, 1 Sep 2023 09:24:22 +0200 Subject: [PATCH 2/4] when showing LossMaps, the user can decide if: * to superimpose loss coordinates to aperture plots; * compute longitudinal histogram of losses (regular LM) or to simply show loss share. In the latter case, the histogram has a color coding assigning to each loss location its own color. --- MADX-tracking/PlotLossMap.m | 49 ++++++++++++++++++++++++++++--------- MADX-tracking/ShowLossMap.m | 38 +++++++++++++++++----------- 2 files changed, 61 insertions(+), 26 deletions(-) diff --git a/MADX-tracking/PlotLossMap.m b/MADX-tracking/PlotLossMap.m index 251e52a..c242e75 100644 --- a/MADX-tracking/PlotLossMap.m +++ b/MADX-tracking/PlotLossMap.m @@ -1,10 +1,10 @@ -function [hh,cc] = PlotLossMap(losses,indices,dS,maxS,minS,Nbins) +function [hh,cc] = PlotLossMap(losses,indices,dS,maxS,minS,lHist,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,indices,dS,maxS,minS,Nbins) +% [hh,cc] = PlotLossMap(losses,indices,dS,maxS,minS,lHist,Nbins) % % input arguments: % losses: table of losses. For the format of the table, please see @@ -15,6 +15,8 @@ % dS: bin width [m]; % maxS: max value of s-coordinate [m]; % minS: min value of s-coordinate [m]; +% lHist: compute longitudinal histogram of losses (boolean); +% otherwise, simply return loss counts; % Nbins: number of bins; % % output arguments: @@ -35,6 +37,7 @@ if ( exist('minS','var') ) minSUsr=minS; end + if (~exist("lHist","var")), lHist=true; end if ( exist('Nbins','var') ) dsUsr=(maxSUsr-minSUsr)/Nbins; end @@ -45,22 +48,44 @@ colS=mapping(find(strcmp(colNames,'s'))); % compute histogram - edges=[minSUsr:dsUsr:maxSUsr]; - if (ismissing(indices)) - [hh,cc] = histcounts(losses{colS},edges); + if (lHist) + % compute longitudinal histogram of losses + edges=[minSUsr:dsUsr:maxSUsr]; + 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 + bar(edges,hh/sum(hh)/dsUsr,1); + ylabel('pdf [m^{-1}]'); + ylim([0.5/(dsUsr*sum(hh)) 2*max(hh)/(sum(hh)*dsUsr)]); else - [hh,cc] = histcounts(losses{colS}(indices),edges); + % compute loss shares + if (ismissing(indices)) + [cc,jj,kk]=unique(losses{mapping(strcmpi(colNames,"S"))}); + else + [cc,jj,kk]=unique(losses{mapping(strcmpi(colNames,"S"))}(indices)); + end + nUniques=length(cc); + hh=accumarray(kk,1); + cm=colormap(parula(nUniques)); + yMin=min(hh)/sum(hh); + % plot histogram + for ii=1:nUniques + if (ii>1), hold on; end + plot([cc(ii) cc(ii)],[yMin/10 hh(ii)/sum(hh)]*100,"Color",cm(ii,:),"LineWidth",1); + end + % bb=bar(cc,hh/sum(hh)*100,1,"EdgeColor","none"); + % bb.CData=cm; + ylabel('[%]'); + ylim([yMin 1]*100); end - edges = cc(2:end) - 0.5*(cc(2)-cc(1)); - - % plot histogram - bar(edges,hh/sum(hh)/dsUsr,1); % additionals xlabel('s [m]'); - ylabel('pdf [m^{-1}]'); xlim([minSUsr maxSUsr]); - ylim([0.5/(dsUsr*sum(hh)) 2*max(hh)/(sum(hh)*dsUsr)]); grid on; end \ No newline at end of file diff --git a/MADX-tracking/ShowLossMap.m b/MADX-tracking/ShowLossMap.m index f24b0ac..8d7c146 100644 --- a/MADX-tracking/ShowLossMap.m +++ b/MADX-tracking/ShowLossMap.m @@ -1,4 +1,4 @@ -function [hh,cc]=ShowLossMap(losses,indices,dS,geometry,allAperH,allAperV) +function [hh,cc]=ShowLossMap(losses,indices,dS,geometry,allAperH,allAperV,lSupImp,lHist) % ShowLossMap embed the longitudinal distribution of losses (histogram) % into a nice plot, possibly with lattice structure and % aperture model @@ -15,6 +15,9 @@ % geometry: (thick) lens lattice; % allAperH: horizontal aperture model; % allAperV: vertical aperture model; +% lSupImp: superimpose loss coordinates to aperture plot (boolean); +% lHist: compute longitudinal histogram of losses (boolean); +% otherwise, simply return loss counts; % % output arguments: % hh: histogram counts; @@ -23,6 +26,9 @@ % See also ReadLosses, GetVariablesAndMappingParticleData, PlotLossMap, % GetColumnsAndMappingTFS, PlotLattice, PlotAperture. + if (~exist("lSupImp","var")), lSupImp=true; end + if (~exist("lHist","var")), lHist=true; end + nPlots=1; if ( exist('geometry','var') ) nPlots=nPlots+1; @@ -71,12 +77,14 @@ tmpAx=subplot(nPlots,1,iPlot); axs=[ axs tmpAx ]; PlotAperture(allAperH{1},allAperH{2},allAperH{3},allAperH{4},allAperH{5}); - hold on; - % losses - if (ismissing(indices)) - plot(losses{colS},losses{colX},'r.'); - else - plot(losses{colS}(indices),losses{colX}(indices),'r.'); + if (lSupImp) + hold on; + % losses + if (ismissing(indices)) + plot(losses{colS},losses{colX},'r.'); + else + plot(losses{colS}(indices),losses{colX}(indices),'r.'); + end end title(sprintf('H plane')); grid on; xlim([minS maxS]); @@ -88,12 +96,14 @@ tmpAx=subplot(nPlots,1,iPlot); axs=[ axs tmpAx ]; PlotAperture(allAperV{1},allAperV{2},allAperV{3},allAperV{4},allAperV{5}); - hold on; - % losses - if (ismissing(indices)) - plot(losses{colS},losses{colY},'r.'); - else - plot(losses{colS}(indices),losses{colY}(indices),'r.'); + if (lSupImp) + hold on; + % losses + if (ismissing(indices)) + plot(losses{colS},losses{colY},'r.'); + else + plot(losses{colS}(indices),losses{colY}(indices),'r.'); + end end title(sprintf('V plane')); grid on; xlim([minS maxS]); @@ -103,7 +113,7 @@ iPlot=iPlot+1; tmpAx=subplot(nPlots,1,iPlot); axs=[ axs tmpAx ]; - [hh,cc] = PlotLossMap(losses,indices,dsUsr,maxS,minS); + [hh,cc] = PlotLossMap(losses,indices,dsUsr,maxS,minS,lHist); set(tmpAx, 'YScale', 'log'); % manual treatment of yticklabels... % set(tmpAx, 'YTick', [min(ylim):10:max(ylim)]); From 7463c14a664d0199bbab60ef2d9d2e416c33bee3 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Fri, 1 Sep 2023 09:26:37 +0200 Subject: [PATCH 3/4] When preparing the scatter plot, it is possible to color initial particles according to properties of the losses --- MADX-tracking/ScatterPlot.m | 89 +++++++++++++++++++++++++++++-------- 1 file changed, 71 insertions(+), 18 deletions(-) diff --git a/MADX-tracking/ScatterPlot.m b/MADX-tracking/ScatterPlot.m index 69c1d9f..f48022c 100644 --- a/MADX-tracking/ScatterPlot.m +++ b/MADX-tracking/ScatterPlot.m @@ -1,48 +1,101 @@ -function ScatterPlot(data,indices,xIDs,yIDs,whichData,currTitle) +function ScatterPlot(data,indices,xNAMEs,yNAMEs,whichData,currTitle,InitialCoords,recogParamName,rPrec) % ScatterPlot Prepare scatter plots % % The user decides what to show % -% ScatterPlot(data,indices,xIDs,yIDs,whichData,currTitle) +% ScatterPlot(data,indices,xNAMEs,yNAMEs,whichData,currTitle) % % input arguments: % data: table with particle coordinates; % indices: index of particles to be plotted (eg to apply some filtering); -% xIDs, yIDs: ID of the quantity to be plotted: +% xNAMEs, yNAMEs: names of the quantity to be plotted: % whichData: format of the table; % currTitle: title of the plot; +% InitialCoords (optional): initial coordinates of particles; +% recogParamName (optional): name of parameter which should be used to +% recognise/label initial conditions; +% rPrec: precision in recognising a given value of recogParamName; +% +% The presence of the InitialCoords, recogParamName and rPrec (~missing()) +% triggers the recognition mode, i.e. points in the scatter plot are +% correlated to unique values of recogParamName; % % See also GetVariablesAndMapping. - if ( length(xIDs)~=length(yIDs) ) + if ( length(xNAMEs)~=length(yNAMEs) ) error('...number of x-axes different from number of y-axes!'); - return end - [nRows,nCols]=GetNrowsNcols(length(xIDs)); + if (~exist("InitialCoords","var")), InitialCoords=missing(); end + if (~exist("recogParamName","var")), recogParamName=missing(); end + if (~exist("rPrec","var")), rPrec=missing(); end + if (all(ismissing(InitialCoords),"all") & all(ismissing(recogParamName),"all") & all(ismissing(rPrec),"all") ) + lRecog=false; + else + lRecog=true; + end + if (lRecog && ismissing(recogParamName)) + recogParamName="S"; + end + if (lRecog && ismissing(rPrec)) + rPrec=1E-3; + end + [nRows,nCols]=GetNrowsNcols(length(xNAMEs)); [ colNames, colUnits, colFacts, mapping ] = GetColumnsAndMappingTFS(whichData); + if (lRecog) + [ colNamesStart, colUnitsStart, colFactsStart, mappingStart ] = GetColumnsAndMappingTFS("starting"); + [uVals,ia,~]=unique(data{mapping(strcmpi(colNames,recogParamName))}); % unique loss locations + uNames=strrep(string(data{mapping(strcmpi(colNames,"ELEMENT"))}(ia)),'"',""); % corresponding names + nUniques=length(uVals); + end actualTitle=sprintf('%s - Scatter plot - %s',whichData,currTitle); ff=figure('Name',actualTitle,'NumberTitle','off'); - for ii=1:length(xIDs) - iColX=mapping(strcmpi(colNames,xIDs(ii))); - iColY=mapping(strcmpi(colNames,yIDs(ii))); + if (lRecog), cm=colormap(parula(nUniques)); end + for ii=1:length(xNAMEs) axt=subplot(nRows,nCols,ii); - if (ismissing(indices)) - if (iscell(data)) - plot(data{iColX}*colFacts(iColX),data{iColY}*colFacts(iColY),'r.'); + if (lRecog) + iColX=mappingStart(strcmpi(colNamesStart,xNAMEs(ii))); + iColY=mappingStart(strcmpi(colNamesStart,yNAMEs(ii))); + % - starting conditions + if (iscell(InitialCoords)) + plot(InitialCoords{iColX}*colFactsStart(iColX),InitialCoords{iColY}*colFactsStart(iColY),'k.'); else - plot(data(:,iColX)*colFacts(iColX),data(:,iColY)*colFacts(iColY),'r.'); + plot(InitialCoords(:,iColX)*colFactsStart(iColX),InitialCoords(:,iColY)*colFactsStart(iColY),'k.'); + end + % - labelling + for jj=1:nUniques + uMin=uVals(jj)-rPrec; uMax=uVals(jj)+rPrec; + indeces=(uMin<=data{mapping(strcmpi(colNames,recogParamName))} & data{mapping(strcmpi(colNames,recogParamName))}<=uMax); + IDs=data{mapping(strcmpi(colNames,"NUMBER"))}(indeces); + hold on; + if (iscell(InitialCoords)) + plot(InitialCoords{iColX}(IDs)*colFactsStart(iColX),InitialCoords{iColY}(IDs)*colFactsStart(iColY),'.',"Color",cm(jj,:)); + else + plot(InitialCoords(IDs,iColX)*colFactsStart(iColX),InitialCoords(IDs,iColY)*colFactsStart(iColY),'.',"Color",cm(jj,:)); + end end + xlabel(sprintf('%s [%s]',colNamesStart(iColX),colUnitsStart(iColX))); + ylabel(sprintf('%s [%s]',colNamesStart(iColY),colUnitsStart(iColY))); else - if (iscell(data)) - plot(data{iColX}(indices)*colFacts(iColX),data{iColY}(indices)*colFacts(iColY),'r.'); + iColX=mapping(strcmpi(colNames,xNAMEs(ii))); + iColY=mapping(strcmpi(colNames,yNAMEs(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 - plot(data(indices,iColX)*colFacts(iColX),data(indices,iColY)*colFacts(iColY),'r.'); + 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))); end - xlabel(sprintf('%s [%s]',colNames(iColX),colUnits(iColX))); - ylabel(sprintf('%s [%s]',colNames(iColY),colUnits(iColY))); grid on; end sgtitle(actualTitle); From 920acd2beb989c1bd04b6d35b736f6edaa0d24ad Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Fri, 1 Sep 2023 09:29:46 +0200 Subject: [PATCH 4/4] actual parsing of multiple loss files --- MADX-tracking/ReadLosses.m | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/MADX-tracking/ReadLosses.m b/MADX-tracking/ReadLosses.m index 4fa6671..ad5e9ae 100644 --- a/MADX-tracking/ReadLosses.m +++ b/MADX-tracking/ReadLosses.m @@ -52,7 +52,7 @@ iCol=mapping(strcmpi(colNames,"TURN")); tmpLosses{iCol}=tmpLosses{iCol}+delays(kk); end - if (~ismissing(nInitial) && ii>1) + if (cumInitial~=0 && ii>1) iCol=mapping(strcmpi(colNames,"NUMBER")); tmpLosses{iCol}=tmpLosses{iCol}+cumInitial(ii-1); end @@ -60,7 +60,9 @@ if (length(losses)==0) losses=tmpLosses; else - losses{nTotLosses+1:nTotLosses+nLosses}=tmpLosses; + for jj=1:size(losses,2) + losses{jj}(nTotLosses+1:nTotLosses+nLosses)=tmpLosses{jj}; + end end nTotLosses=nTotLosses+nLosses; end