Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Changes necessary to perform a Beam tomography #65

Open
wants to merge 4 commits into
base: main
Choose a base branch
from
Open

[WIP] Changes necessary to perform a Beam tomography #65

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
2e34400
`GetColumnsAndMappingTFS.m` now also handles start contiditions of PT…
Aug 31, 2023
1614282
when showing LossMaps, the user can decide if:
Sep 1, 2023
7463c14
When preparing the scatter plot, it is possible to color initial part…
Sep 1, 2023
920acd2
actual parsing of multiple loss files
Sep 1, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
12 changes: 8 additions & 4 deletions MADX-optics/GetColumnsAndMappingTFS.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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" ...
Expand All @@ -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
Expand All @@ -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" ];
Expand All @@ -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" ...
Expand All @@ -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?');
Expand Down
49 changes: 37 additions & 12 deletions MADX-tracking/PlotLossMap.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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:
Expand All @@ -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
Expand All @@ -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
6 changes: 4 additions & 2 deletions MADX-tracking/ReadLosses.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -52,15 +52,17 @@
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
nLosses=length(tmpLosses{1});
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
Expand Down
89 changes: 71 additions & 18 deletions MADX-tracking/ScatterPlot.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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);
Expand Down
38 changes: 24 additions & 14 deletions MADX-tracking/ShowLossMap.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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;
Expand All @@ -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;
Expand Down Expand Up @@ -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]);
Expand All @@ -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]);
Expand All @@ -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)]);
Expand Down