changes necessary for performing a Beam tomography

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

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;

MADX-optics/ParseTfsTable.m

@@ -23,34 +23,52 @@
 %   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; % 52; % 48;
+            case "PTC_TWISS"
+                nHeader=90;
+            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);

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=[];

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