From 3df1dc2b00ecea09dc1c6ab043cdfe8990da4f1a Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Mon, 21 Nov 2022 17:19:15 +0100 Subject: [PATCH 01/23] adding a couple of useful functions: * one for determining, adding and checking the path to Reference Machine Parameters repo * one for decoding and checking machine configuration --- operations/DecodeConfig.m | 102 +++++++++++++++++++++++++++++++++ operations/PathToMachineRefs.m | 37 ++++++++++++ 2 files changed, 139 insertions(+) create mode 100644 operations/DecodeConfig.m create mode 100644 operations/PathToMachineRefs.m diff --git a/operations/DecodeConfig.m b/operations/DecodeConfig.m new file mode 100644 index 0000000..f851125 --- /dev/null +++ b/operations/DecodeConfig.m @@ -0,0 +1,102 @@ +function [machine,beamPart,focus,config]=DecodeConfig(myConfig) +% DecodeConfig(myConfig) decodes the configuration in its parts, after +% some checking; +% input (all optional): +% - myConfig (string): machine configuration to be found. It is a string +% made of substrings separated by comas, no empty spaces. +% The string should look like: ",,,"; +% e.g. "SYNCHRO,PROTON,TM" +% The user can specify a subset of values; for the others, the default +% values will be used; +% +% ouput: +% - machine (string): beam line (default: SYNCHRO); +% - beamPart (string): transported particle (default: PROTON); +% - focus (string): dimension of focus (default: FG); +% - config (string): extraction method (default: TM); +% +% See also ReturnDefFile. + + stripMeAll=[" " "-" "_"]; + + % - default values (focus depends on the particle) + if (~exist("myConfig","var")), myConfig=missing(); end + machine="SYNCHRO"; beamPart="PROTON"; config="TM"; + % - actual requests by user + if (~ismissing(myConfig) && strlength(myConfig)>0) + myInfo=split(myConfig,","); + if (length(myInfo)>=1) + if ( strlength(myInfo(1))>0 ) + machine=upper(myInfo(1)); % capitalize + end + end + if (length(myInfo)>=2) + if ( strlength(myInfo(2))>0 ) + beamPart=upper(myInfo(2)); % capitalize + end + end + if (length(myInfo)>=3) + if ( strlength(myInfo(3))>0 ) + config=upper(myInfo(3)); % capitalize + end + end + if (length(myInfo)>=4) + if ( strlength(myInfo(4))>0 ) + focus=upper(myInfo(4)); % capitalize + end + end + end + % - take into account nick-names, typos, etc... + switch erase(machine,stripMeAll) % erase characters + case {"LINEZ","SALA1"} + machine="Sala1"; + case {"LINEV","SALA2V"} + machine="Sala2V"; + case {"LINEU","SALA2H"} + machine="Sala2H"; + case {"LINET","SALA3"} + machine="Sala3"; + case {"ISO1","XPRX1","X1"} + machine="Iso1"; + case {"ISO2","XPRX2","X2"} + machine="Iso2"; + case {"ISO3","XPRX3","X3"} + machine="Iso3"; + case {"ISO4","XPRX4","X4"} + machine="Iso4"; + case {"SINCRO","SYNCHRO"} + machine="SYNCHRO"; + otherwise + error("machine %s NOT recognised!",machine); + end + switch erase(beamPart,stripMeAll) % erase characters + case {"PROTON","PROTONI","PROT","P"} + beamPart="P"; + if (~exist("focus","var")), focus="FG"; end + case {"CARBON","CARBONIO","CARB","C"} + beamPart="C"; + if (~exist("focus","var")), focus="FP"; end + otherwise + error("particle %s NOT recognised!",beamPart); + end + switch erase(config,stripMeAll) % erase characters + case {"BETATRONE","BETATRON","BETA","BET"} + config="Betatrone"; + case {"RFKO","RFKNOCKOUT"} + config="RFKO"; + case "TM" + config="Betatrone"; % TM=Betatrone + otherwise + error("configuration %s NOT recognised!",config); + end + switch erase(focus,stripMeAll) % erase characters + case {"FG","GRANDE"} + focus="FG"; + case {"FP","PICCOLO"} + focus="FP"; + otherwise + error("focus %s NOT recognised!",focus); + end + if (strcmp(machine,"SYNCHRO")), focus=""; end + +end diff --git a/operations/PathToMachineRefs.m b/operations/PathToMachineRefs.m new file mode 100644 index 0000000..ee4902f --- /dev/null +++ b/operations/PathToMachineRefs.m @@ -0,0 +1,37 @@ +function MRpath=PathToMachineRefs() +% PathToMachineRefs returns the path to the MachineRefs clone. +% This is searched for in the same path where +% MatLabTools is located. +% If the path is not found in the path env var, +% then it is automatically added. + MRrepoName="MachineRefs"; + MLTrepoName="MatLabTools"; + if (contains(path,MRrepoName)) + allPaths=string(split(path,";")); + MRpaths=allPaths(contains(allPaths,MRrepoName)); + if (isempty(MRpaths)) + error("Something wrong with getting path to %s",MRrepoName); + end + MRpath=extractBetween(MRpaths(1),1,MRrepoName,"Boundaries","inclusive"); + if (~isfolder(MRpath)) + error("Cannot access path to %s: %s",MRrepoName,MRpath); + end + else + warning("path to %s not loaded. Building it...",MRrepoName); + if (~contains(path,MLTrepoName)) + error("Path to %s NOT loaded",MLTrepoName); + else + allPaths=string(split(path,";")); + MLTpaths=allPaths(contains(allPaths,MLTrepoName)); + if (isempty(MLTpaths)) + error("Something wrong with getting path to %s",MLTrepoName); + end + MLTpath=extractBetween(MLTpaths(1),1,MLTrepoName,"Boundaries","inclusive"); + if (~isfolder(MLTpath)) + error("Cannot access path to %s: %s",MLTrepoName,MLTpath); + end + MRpath=strcat(extractBetween(MLTpaths(1),1,MLTrepoName),"/",MRrepoName); + addpath(genpath(MRpath)); + end + end +end From 70ef913d6ad04471bc280e5020039c7c4448610f Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Tue, 22 Nov 2022 11:40:29 +0100 Subject: [PATCH 02/23] bug fixes in decoding configuration --- operations/DecodeConfig.m | 31 +++++++++++++++---------------- 1 file changed, 15 insertions(+), 16 deletions(-) diff --git a/operations/DecodeConfig.m b/operations/DecodeConfig.m index f851125..9cebff2 100644 --- a/operations/DecodeConfig.m +++ b/operations/DecodeConfig.m @@ -21,31 +21,28 @@ % - default values (focus depends on the particle) if (~exist("myConfig","var")), myConfig=missing(); end - machine="SYNCHRO"; beamPart="PROTON"; config="TM"; + machine="SYNCHRO"; beamPart="PROTON"; config="TM"; focus=""; % - actual requests by user if (~ismissing(myConfig) && strlength(myConfig)>0) myInfo=split(myConfig,","); if (length(myInfo)>=1) - if ( strlength(myInfo(1))>0 ) - machine=upper(myInfo(1)); % capitalize - end + if ( strlength(myInfo(1))>0 ), machine=myInfo(1); end end if (length(myInfo)>=2) - if ( strlength(myInfo(2))>0 ) - beamPart=upper(myInfo(2)); % capitalize - end + if ( strlength(myInfo(2))>0 ), beamPart=myInfo(2); end end if (length(myInfo)>=3) - if ( strlength(myInfo(3))>0 ) - config=upper(myInfo(3)); % capitalize - end + if ( strlength(myInfo(3))>0 ), config=myInfo(3); end end if (length(myInfo)>=4) - if ( strlength(myInfo(4))>0 ) - focus=upper(myInfo(4)); % capitalize - end + if ( strlength(myInfo(4))>0 ), focus=myInfo(4); end end end + % - capitalize info + machine=upper(machine); + beamPart=upper(beamPart); + config=upper(config); + focus=upper(focus); % - take into account nick-names, typos, etc... switch erase(machine,stripMeAll) % erase characters case {"LINEZ","SALA1"} @@ -72,10 +69,10 @@ switch erase(beamPart,stripMeAll) % erase characters case {"PROTON","PROTONI","PROT","P"} beamPart="P"; - if (~exist("focus","var")), focus="FG"; end + if (~exist("focus","var") || strlength(focus)==0), focus="FG"; end case {"CARBON","CARBONIO","CARB","C"} beamPart="C"; - if (~exist("focus","var")), focus="FP"; end + if (~exist("focus","var") || strlength(focus)==0), focus="FP"; end otherwise error("particle %s NOT recognised!",beamPart); end @@ -95,7 +92,9 @@ case {"FP","PICCOLO"} focus="FP"; otherwise - error("focus %s NOT recognised!",focus); + if (strlength(focus)>0) + error("focus %s NOT recognised!",focus); + end end if (strcmp(machine,"SYNCHRO")), focus=""; end From 3a7eb76e29533537b543118ce39fef31f115ed50 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Tue, 22 Nov 2022 11:49:11 +0100 Subject: [PATCH 03/23] AcquireLGENValues.m now uses the MachineRef repo --- CreateQAfiles.m | 13 ++-- operations/AcquireLGENValues.m | 129 +++++---------------------------- 2 files changed, 25 insertions(+), 117 deletions(-) diff --git a/CreateQAfiles.m b/CreateQAfiles.m index 8715e1d..cb1d173 100644 --- a/CreateQAfiles.m +++ b/CreateQAfiles.m @@ -10,29 +10,32 @@ %% include libraries % - include Matlab libraries -pathToLibrary=".\"; +pathToLibrary="./"; +addpath(genpath(pathToLibrary)); +pathToLibrary="../MachineRefs"; addpath(genpath(pathToLibrary)); %% settings % ------------------------------------------------------------------------- % USER's input data -kPath="S:\Accelerating-System\Accelerator-data"; +kPath="P:\Accelerating-System\Accelerator-data"; % kPath="K:"; -beamPart="CARBON"; machine="ISO3"; +beamPart="CARBON"; config="TM"; % select configuration: TM, RFKO % ------------------------------------------------------------------------- %% parse DBs +myConfig=sprintf("%s,%s,%s",machine,beamPart,config); % - get PS mapping -clear PSmapping; PSmapping=readtable("PSmapping.xlsx"); +clear PSmapping; FullFileName=ReturnDefFile("PSmapping",myConfig); PSmapping=readtable(FullFileName); % - get TM values clear cyCodesTM rangesTM EksTM BrhosTM currentsTM fieldsTM kicksTM psNamesTM FileNameCurrentsTM magNamesTM ; -[cyCodesTM,rangesTM,EksTM,BrhosTM,currentsTM,fieldsTM,kicksTM,psNamesTM,FileNameCurrentsTM]=AcquireLGENValues(beamPart,machine,config); +[cyCodesTM,rangesTM,EksTM,BrhosTM,currentsTM,fieldsTM,kicksTM,psNamesTM,FileNameCurrentsTM]=AcquireLGENValues(myConfig); psNamesTM=string(psNamesTM); cyCodesTM=upper(string(cyCodesTM)); magNamesTM=MagNames2LGENnames(psNamesTM,true,PSmapping); diff --git a/operations/AcquireLGENValues.m b/operations/AcquireLGENValues.m index 1df8d2b..8aa7e46 100644 --- a/operations/AcquireLGENValues.m +++ b/operations/AcquireLGENValues.m @@ -1,10 +1,19 @@ %% main function -function [cyCodes,ranges,Eks,Brhos,currents,fields,kicks,psNames,FileNameCurrents]=AcquireLGENValues(beamPartIn,machineIn,configIn,filters,LGENsCheck) +function [cyCodes,ranges,Eks,Brhos,currents,fields,kicks,psNames,FileNameCurrents]=AcquireLGENValues(myConfig,path2Files,filters,LGENsCheck) fprintf("Getting current values...\n"); - % filter PSs where all CyCo show "NaN" or "0" + %% input arguments + [machine,beamPart,focus,config]=DecodeConfig(myConfig); + fprintf("...for: machine=%s; beamPart=%s; config=%s;\n",machine,beamPart,config); + if (~exist("path2Files","var")), path2Files=missing(); end + if (ismissing(path2Files)) + path2Files=ReturnDefFile("LGEN",myConfig); + end + if (length(path2Files)==1) + path2Files(2)=ReturnDefFile("BRHO",myConfig); + end if ( ~exist('filters','var') ), filters=["NaN" "0"]; end if ( exist('LGENsCheck','var') ) doVisualCheck=true; @@ -13,32 +22,7 @@ LGENsCheck=missing; end - % processing - beamPart=upper(beamPartIn); - machine=upper(machineIn); - config=upper(configIn); - - % preliminary checks - if ( ~strcmp(config,"TM") & ~strcmp(config,"RFKO") ) - error("unrecognised config: %s - available only TM and RFKO!",config); - end - if ( ~strcmp(beamPart,"PROTON") & ~strcmp(beamPart,"CARBON") ) - error("unrecognised beam particle: %s - available only PROTON and CARBON!",beamPart); - end - if ( ~strcmp(machine,"SYNCHRO") & ~strcmp(machine,"LINET") & ~strcmp(machine,"SALA3") ... - & ~strcmp(machine,"LINEU") & ~strcmp(machine,"SALA2H") ... - & ~strcmp(machine,"LINEV") & ~strcmp(machine,"SALA2V") ... - & ~strcmp(machine,"LINEZ") & ~strcmp(machine,"SALA1") ... - & ~strcmp(machine,"XPRX1") & ~strcmp(machine,"ISO1") ... - & ~strcmp(machine,"XPRX2") & ~strcmp(machine,"ISO2") ... - & ~strcmp(machine,"XPRX3") & ~strcmp(machine,"ISO3") ... - & ~strcmp(machine,"XPRX4") & ~strcmp(machine,"ISO4") ... - ) - error("unrecognised machine: %s - available only SYNCHRO, LINEZ/SALA1, LINEV/SALA2V, LINEU/SALA2H, LINET/SALA3, XPRX1/ISO1, XPRX2/ISO2, XPRX3/ISO3 and XPRX4/ISO4!",machine); - end - fprintf("...for: machine=%s; beamPart=%s; config=%s;\n",machine,beamPart,config); - - % do the job + %% do the job switch machine case "SYNCHRO" switch beamPart @@ -81,91 +65,12 @@ CyCoData={CyCoData{:,1} ; CyCoData{:,2} ; CyCoData{:,4} ; temp{:,1} }'; buffer = vertcat( CyCoData{:,1} ) ; % extract only first four digits of cyco CyCoData(:,1) = cellstr(buffer(:,1:4)) ; % - case {"LINEZ","SALA1","LINEV","SALA2V","LINEU","SALA2H","LINET","SALA3",... - "XPRX1","ISO1","XPRX2","ISO2","XPRX3","ISO3","XPRX4","ISO4"} - switch beamPart - case "PROTON" - % BUILD TABLEs WITH CyCo, Range, Energy and Brho - FileName="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\MeVvsCyCo_P.xlsx"; - CyCoData = GetOPDataFromTables(FileName,"Sheet1"); - FileName="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\EvsBro_P.xlsx"; - BrhoData = GetOPDataFromTables(FileName,"Sheet1"); - % PARSE FILE WITH CURRENTS AT FT - columns in final cell array: - switch machine - case {"LINEZ","SALA1"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\HEBT\Protoni\ProtoniSala1\Protoni_Sala1_2022-03-09.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents,"09.03.2022 - 10.27"); - case {"LINEV","SALA2V"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\HEBT\Protoni\ProtoniSala2V\Protoni_Sala2V_2021-02-13.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents,"27.01.2022 - 09.08"); - case {"LINEU","SALA2H"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\HEBT\Protoni\ProtoniSala2H\Protoni_Sala2H_2021-08-09.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents,"09.08.2021 - 14.51"); - case {"LINET","SALA3"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\HEBT\Protoni\ProtoniSala3\Protoni_Sala3_2021-08-11.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents,"11.08.2021 - 08.59"); - case {"XPRX1","ISO1"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\XPR\XPR1\Protoni\FuocoGrande\LGEN_X1_P_FG_10Luglio2019.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents); - case {"XPRX2","ISO2"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\XPR\XPR2\Protoni\FuocoGrande\LGEN_X2_P_FG_10Luglio2019.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents); - case {"XPRX3","ISO3"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\XPR\XPR3\Protoni\FuocoGrande\LGEN_X3_P_FG_10Luglio2019.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents); - case {"XPRX4","ISO4"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\XPR\XPR4\Protoni\FuocoGrande\LGEN_X4_P_FG_10Luglio2019.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents); - end - case "CARBON" - % BUILD TABLEs WITH CyCo, Range, Energy and Brho - FileName="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\MeVvsCyCo_C.xlsx"; - CyCoData = GetOPDataFromTables(FileName,"Sheet1"); - FileName="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\EvsBro_C.xlsx"; - BrhoData = GetOPDataFromTables(FileName,"Sheet1"); - % PARSE FILE WITH CURRENTS AT FT - columns in final cell array: - switch machine - case {"LINEZ","SALA1"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\HEBT\Carbonio\lineaZ\fuocopiccolo\Carbonio_Sala1_FromRepoNovembre2020.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents,"09.11.2020 - 10.10"); - case {"LINEV","SALA2V"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\HEBT\Carbonio\lineaV\FuocoPiccolo\Carbonio_Sala2V_FromRepo.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents,"21.08.2019 - 12.11"); - case {"LINEU","SALA2H"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\HEBT\Carbonio\lineaU\FuocoPiccolo\Carbonio_Sala2H_FromRepo.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents,"21.08.2019 - 12.04"); - case {"LINET","SALA3"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\HEBT\Carbonio\lineaT\fuocopiccolo\Carbonio_Sala3_FromRepoNovembre2020.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents,"09.11.2020 - 10.11"); - case {"XPRX1","ISO1"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\XPR\XPR1\Carbonio\FuocoPiccolo\LGEN_X1_C_FP.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents); - case {"XPRX2","ISO2"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\XPR\XPR2\Carbonio\FuocoPiccolo\LGEN_X2_C_FP.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents); - case {"XPRX3","ISO3"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\XPR\XPR3\Carbonio\FuocoPiccolo\LGEN_X3_C_FP_X2fit.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents); - case {"XPRX4","ISO4"} - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\XPR\XPR4\Carbonio\FuocoPiccolo\LGEN_X4_C_FP.xlsx"; - currentData = GetOPDataFromTables(FileNameCurrents); - end - otherwise - error("no source of data available for %s %s %s",machine,beamPart,config); - end % switch: LGEN, LINET/SALA3/LINEU/SALA2H/LINEV/SALA2V/LINEZ/SALA1, beamPart - % - get CyCo (col 2 []), range (col 3 [mm]) and Energy (col 1 [MeV/n]) - columns in cell array - CyCoData = CyCoData(2:end,:); % remove header - % - get Brho (col 2 [Tm]) and Energy (col 1 [MeV/n]) - columns in cell array: - BrhoData = BrhoData(2:end,:); % remove header - % - get common ranges - [commonRanges,iCRa,iCRb]=intersect(cell2mat(CyCoData(:,3)),cell2mat(BrhoData(:,1))); - CyCoData=CyCoData(iCRa,:); - BrhoData=BrhoData(iCRb,:); - % - make a unique table: 1=CyCo[], 2=range[mm], 3=Energy[MeV/n], 4=Brho[Tm] - columns in final cell array - CyCoData={CyCoData{:,2} ; CyCoData{:,3} ; CyCoData{:,1} ; BrhoData{:,2} }'; otherwise - error("no source of data available for %s %s %s",machine,beamPart,config); - end % switch: machine + currentData = GetOPDataFromTables(path2Files(1)); FileNameCurrents=path2Files(1); + CyCoData = GetOPDataFromTables(path2Files(2)); + % - make a unique table: 1=CyCo[], 2=range[mm], 3=Energy[MeV/n], 4=Brho[Tm] - columns in final cell array + CyCoData=CyCoData(2:end,1:4); + end % get currents and PSs to be crunched % - nRows: number of power supplies + a header From 1e378b9ebdc017e11e50ff838c6f80a5276b3555 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Tue, 22 Nov 2022 14:16:20 +0100 Subject: [PATCH 04/23] AcquireLGENValues.m and AcquireLGENValues.m now crunch also SYNCHRO PSs --- operations/AcquireLGENValues.m | 52 +++------------------------------ operations/MagNames2LGENnames.m | 22 ++++++++++---- 2 files changed, 20 insertions(+), 54 deletions(-) diff --git a/operations/AcquireLGENValues.m b/operations/AcquireLGENValues.m index 8aa7e46..b6efde2 100644 --- a/operations/AcquireLGENValues.m +++ b/operations/AcquireLGENValues.m @@ -23,54 +23,10 @@ end %% do the job - switch machine - case "SYNCHRO" - switch beamPart - case "PROTON" - % BUILD TABLE WITH CyCo, Range, Energy and Brho - % - get CyCo (col 1 []), range (col 2 [mm]) and Energy (col 4 [MeV/n]) - columns in cell array - FileName="P:\Accelerating-System\Accelerator-data\Area dati MD\00Rampe\MatlabRampGen2.8\INPUT\CSV-TRATTAMENTI\Protoni.csv"; - CyCoData = GetOPDataFromTables(FileName); - % - build array of values of Brho (as in RampGen) - mp = 938.255; An = 1; Zn = 1; - % PARSE FILE WITH CURRENTS AT FT - columns in final cell array: - % - nRows: number of power supplies + a header - % - nColumns: number of cycle codes + 2 (PS name + property) - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\Sincro\CorrentiFlatTop\ProtoniSincro_2021-02-13.xlsx"; % AMereghetti, 2021-11-19: no longer there! - currentData = GetOPDataFromTables(FileNameCurrents,"Foglio1"); - case "CARBON" - if ( strcmp(config,"RFKO") ) - error("no source of data available for %s %s %s",machine,beamPart,config); - end - % BUILD TABLE WITH CyCo, Range, Energy and Brho - % - get CyCo (col 1 []), range (col 2 [mm]) and Energy (col 4 [MeV/n]) - columns in cell array - FileName="P:\Accelerating-System\Accelerator-data\Area dati MD\00Rampe\MatlabRampGen2.8\INPUT\CSV-TRATTAMENTI\Carbonio.csv"; - CyCoData = GetOPDataFromTables(FileName); - % - build array of values of Brho (as in RampGen) - mp = 931.2225; An = 12; Zn = 6; - % PARSE FILE WITH CURRENTS AT FT - columns in final cell array: - % - nRows: number of power supplies + a header - % - nColumns: number of cycle codes + 2 (PS name + property) - FileNameCurrents="P:\Accelerating-System\Accelerator-data\Area dati MD\00Setting\Sincro\CorrentiFlatTop\CarbonioSincro_2021-02-05.xlsx"; % AMereghetti, 2021-11-19: no longer there! - currentData = GetOPDataFromTables(FileNameCurrents,"Foglio1"); - otherwise - error("no source of data available for %s %s %s",machine,beamPart,config); - end % switch: LGEN, SYNCHRO, beamPart - % continue crunching CyCo data - CyCoData = CyCoData(2:end,:); % remove header - c = 2.99792458e8; % velocità della luce [m/s] - BRO = @(x)(An/Zn)*((mp*sqrt((1 + x/mp).^2 - 1))/c)*10^6; - temp=num2cell(BRO(cell2mat(CyCoData(:,4)))); - % - make a unique table: 1=CyCo[], 2=range[mm], 3=Energy[MeV/n], 4=Brho[Tm] - columns in final cell array - CyCoData={CyCoData{:,1} ; CyCoData{:,2} ; CyCoData{:,4} ; temp{:,1} }'; - buffer = vertcat( CyCoData{:,1} ) ; % extract only first four digits of cyco - CyCoData(:,1) = cellstr(buffer(:,1:4)) ; % - otherwise - currentData = GetOPDataFromTables(path2Files(1)); FileNameCurrents=path2Files(1); - CyCoData = GetOPDataFromTables(path2Files(2)); - % - make a unique table: 1=CyCo[], 2=range[mm], 3=Energy[MeV/n], 4=Brho[Tm] - columns in final cell array - CyCoData=CyCoData(2:end,1:4); - end + currentData = GetOPDataFromTables(path2Files(1)); FileNameCurrents=path2Files(1); + CyCoData = GetOPDataFromTables(path2Files(2)); + % - make a unique table: 1=CyCo[], 2=range[mm], 3=Energy[MeV/n], 4=Brho[Tm] - columns in final cell array + CyCoData=CyCoData(2:end,1:4); % get currents and PSs to be crunched % - nRows: number of power supplies + a header diff --git a/operations/MagNames2LGENnames.m b/operations/MagNames2LGENnames.m index 46d5212..354ebbe 100644 --- a/operations/MagNames2LGENnames.m +++ b/operations/MagNames2LGENnames.m @@ -1,24 +1,34 @@ -function outNames=MagNames2LGENnames(inNames,lInv,PSmapping) +function outNames=MagNames2LGENnames(inNames,lInv,PSmapping,separator) if (~exist("lInv","var")), lInv=false; end % by default, convert Magnet Name into LGEN + if (~exist("separator","var")), separator=","; end outNames=strings(length(inNames),1); if (lInv) % LGEN name --> Magnet name arrayIN=PSmapping.LGEN; arrayOUT=PSmapping.Magnet_NAME; - whatIN="Magnet Name"; + whatTO="Magnet Name"; else % Magnet name --> LGEN name arrayIN=PSmapping.Magnet_NAME; arrayOUT=PSmapping.LGEN; - whatIN="LGEN Name"; + whatTO="LGEN Name"; end % query mapping for ii=1:length(inNames) - jj=find(strcmpi(string(arrayIN),inNames(ii))); - if ( isempty(jj) ), error("%s %s not found in LGEN mapping table!",whatIN,inNames(ii)); end - outNames(ii)=string(arrayOUT(jj)); + if (any(contains(string(arrayIN),separator))) + jj=find(contains(string(arrayIN),inNames(ii))); + else + jj=find(strcmpi(string(arrayIN),inNames(ii))); + end + if ( isempty(jj) ), error("%s %s not found in LGEN mapping table!",whatTO,inNames(ii)); end + if (any(contains(string(arrayOUT),separator))) + myStrings=split(arrayOUT(jj),separator); + outNames(ii)=myStrings(1); + else + outNames(ii)=string(arrayOUT(jj)); + end end end From 73135d8fedeb81bb670da4144820ba6f50e7c86a Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Tue, 22 Nov 2022 14:21:30 +0100 Subject: [PATCH 05/23] removing available mapping. It has been moved to MachineRefs repo --- MeVvsCyCo_P.xlsx | Bin 15282 -> 0 bytes PSmapping.xlsx | Bin 12437 -> 0 bytes 2 files changed, 0 insertions(+), 0 deletions(-) delete mode 100644 MeVvsCyCo_P.xlsx delete mode 100644 PSmapping.xlsx diff --git a/MeVvsCyCo_P.xlsx b/MeVvsCyCo_P.xlsx deleted file mode 100644 index 7fa3748c2c48bdcf55b411a396a4e7e82b66e52e..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 15282 zcmeHuWpErzmbF;2n31%st_sw|TMtuKw ztE0QC^7cKIdGkU>*GVN=a0qk|C=eJB5D*d&l30oNuAm?wv5+7jXdp0Px}uH_Zsrbd zhHBnU=C1mTUiNmx#SmcBg&<(S-~Zq7zt{o;YNHN)%*dUZM}lHq8d1Tb>go`P`tb~? zKf&erNHC-9y`f+O*Q2!6)agZtZ<47S_U&)w0$Xj4%-(^SVGh3nOV4mvpnE)>-A8}-d$x+H>f#J1#iY&9cq5V4*Ukp} zOoxdeYIy)o9ftwJ;l3AYJj4t|${8U-xC5`WZC!OS_1_fI_n-p`Hlp#nV&+w=s_KsS4P0|!z1w^(gbXD0gzl*j=?5e^uu zhA!rIu1t)-KL0O9{}%`5zYM)LSzfW989w4%>N#@cZf-RZSwz-TSgM0WEg(pG1*s{n zgdA_Rn+hLU4L1~AGO#n?d2D%=FYaiBA zBC4k!BkEhY3Xf7Ia{?FA%MTF*Kk&$%%w%DVxSCileexTzBe{9NRadp*v934CcX?0h zZD?lKb1jzBjrQWlppZMFLduSM&pIhNOi_3fsPmcq@@t;=0N1F8o>;^=7mtEBx`%zy&x8W81w_wglp+-86o zKI}ZaD|EWkIr&v-?7c5dcxUr6G;;%Lq`aMJm7w=icIp{teQCy22|4>U%QuTBo}XI; zWbSJ;)Bt>l_8&PDO|S-cikf%CADYEg4S3lz#`3=Ipr;+oQhenM&&*WP=`i|ytp%(jlsQ%3gO!yO2KAANANUuV~kOXv!~ zo4Uk-8;ycXDU4?TgjhctU}}RwdgmAErEK;OgKLWkIOA)QSID+Os6MlOGhuD9?Ms|P z8HXEI6L??7rprX(+DHqZ;D7YkK*v5VcJ6LGLtgqLdWFKLz1CZZ;lir5oa%gl9%8xMPsmuXun6{+@ynFXqM?KN$CX5 zp_Bq)Kuss&-QAXByQ1SJ>#Mr+oWV)>5IRGvFF{Qq;czQ?9Kvm7_oT;Y$r{aNRbu8e z-DDk1&57YX>!@81$kkc1*Dc$0$K^#?d&8i#oKCeunXQwcu(D|vnL}{mP)iM46Z|Pt z(o#!npw#@NHK3ow$pp#LnF6nNv02tn&jf+9^q(0Z2tcW{ z2%JL;Xb=#Df6Sw+o41|0>#xWkQ2#xPomdBg$Y1DmOxl`BXg+Miz$9XVCWCd(pqis? z7OSYyR*ZCeQjayyps%VFVkNQR3v}TK#3wf@IQ1lh z_hiirlu4^dy@b{k%&DuIR%0rIe*bj^x|F^ul^uarjSk^EY~ks<7~iU{krhVmikPwF zJrsUq#C+r@Z58TI?5s*S#Ph!7DBs2PkSFG|m+GTohCC+-v?UF{kti@rjlcK*YKo1! zIKVAfA59PHWcE&NSc$=vd5XDvdCuueSd1I9{|ZwvSN2ICHKF~*6L_Q9N4Ix$36Z;`%XsGvFH1rCpEf@8}&FEZ=KaHV> z%Gsmbs!EwPYe}m2M`NF)et3rxFC{5XA_s6?etA6fTl56x@6JJJY;yS&* z=y_l6n6&BaUK7JI=1Z(zVsti~ZxZJ^0j#*)l=#VcOTB!5mU)?RKl+=3{xChew!+OF zHs08+5I&GEH`^O8uV}~^HvXKgke=a8n9w>LILr_K%A1od|zz-!Mc+Y5Ak5nwta`zw{XTA7=>x&E^r{=Q)O?`jy9&>P*y zj4pl${wnI_n=&5`FRIE#+N!<}Zu0m8YB?kE4*kWq5#PLLch_lhmow*yuL{osF^T?r zJ6QORDVA0Jh3^lBZs7**qZRab419^kvy5HY;0g*dQNxK^7;p{Y%H%KD{Q(0@4e{@) z5(V~_q@%oDYEi%hF|{_AHTEl;Ww6i2I;pS<0iInh(sTXVxknO4QJNr9NFq&_xo~AfqV6 ziqv@;M5Pg{yWvw2jR6<`N-DzA;&JEOzr(iLo3N{5l9_7xm<9!7WG>xl*%0 zXa-M!N%E(%~Q6rTO1x6IS63f7#Nn0a}Q|=w+k{i58_)Qo* z+5st&opO3e%r7!HUQe1InXRENJjEqGlmK6iLnPij=bhUJZ8k7@@PO}RNNQo&>9B71V!H*Z`}RpIvDU$^g7>FVrG2dO12?iETsU<5k$0q) z54yn}k*1@om(_Hp)5MW<7<4l1y>u$msd9Z8|D48?gZELRW{ujuK?GqWpmuz7fwrb; z+#JbK2-C2je?ceb^tL*Zi#WzGuivwdUBgyjEl5S;I+W)#qE(BoOo2#5u08aoEr?-F z>IDugz&Kl9XE1T!M){NZy#SUzh!t^v^qy249P0qfjx(O6zp!%{Y_97&@U&f-764e_ zjpMjqwe^ixo{?ofLG4f#CDu{dCE8 zq@a%W+JJp7)!9JCdcf`e0#S)ze!%tBC6@6Qx6AI=i%}-QKtd=H|MB^^t*=zw{JvLP zHHv|6XU-s9G*0W`cUBVQ?c(uMgp)@0Lo8jprA!lCB^f41<2^JZh#$@>zk{#gnqu}f zeK`Np0ai%?@kOBhe5t)4I;j|RoNMMA_?gp$jvfSoljyMv-wNKYv{oB-TPuGol4$!f z_npWx4n|K^BH+5qv0{MDHDgn6Gi*nU)=bIX)}0$mxT{KBJ8@rp^r}&(9jSh1 zIBkNM4Kc@3kVy(K+m&p@DiM1tnDqu!2VBhc+RP`scDWByMF+xCd-C7fVpg`GPth&K zr!2B77;nM7L1?0<88u1n5%1{40wcLM{ES=ARkp&P+(=tR(I!i^Lk_6OP(JLs_r3Vw zHyNTidCXkU%IXcnr4f!s){XpPT7jkLOLR~}>4^VQ)F*BKfde87 z2b=<|nYnK$fAYh$Q(Rop>fu@~Z?cEs;8=ApyFQcC>)eD#NR@jDjbR_UPYzewgs_jM zN#jDAluqVV42(?APnoQXENh4CRM#*laYd~yi&?uY=8H=&oJP5rwdv;4OG=JcDjR;1 zq#hDy!|h|i6R?03(?8NcYlLv9F`=;gQK0jUQr~jgE^HP6%SAPR?h1G7Rfw!&Z=Chx zJ9EMqx=~vjI`b{lLhnbcCXk#%N29?pT=M>P3&{IvBnBKr>Dm<|Ua~!n0T3D5T{7{3 zfcM-Y{+y`3*nM@X81rm}sAr<^WCu)6O%4qhM_sWqV9gR!hIz{R?Ca*{EsjRl8C zrw-|ZF4#oy66Y(qe?2fH7p9@QQ0S zt#s5%X~44L>}Jxz*?JV<2Q@RTdmj?ic)cSTqrE z^eyZFUp85KusjjBK|JKu6((!F%gVQrfv#IxoNSxZ8OqkZbb0xvYY)00=hRV4Ol3lA zVqutmfY3C27o>~w<*EJ@Bc!&=3=gUF-OP+B5iS+4Rf=>B06}(X)xZl6l55Ajf_Dy1 z{gE>7j3p)J)3-Tu-$7{g2Cv-&HJCU*8eG z&BY%n@s0T4r8(OUAPB#VrKv*?5CVf(^5aX`ZD6P^RW5Dl%rPtw1?bplG_cd^Z+v3N zu!*nM<`V}XKE15Jed*wTY6$de=UEBl_S6pyXzuP7@at-Od)$2Pe!X7k9^bcA9A16d zUp?Inc-wrgdzpP$7ks@I?0oa~dfq-<4hZB6jeM#=+Fy?!j`YOBQe0TS_3WrPyQuFT z?)Ky1TRV-vKi*!5sfZkYIdgq2>R5(bzuO*MJ==Tz;^lRBk4Y5CWa#ket{eYT_xkC1 z!r%E(+uj?}%FTY?%N5I|k!c+hD z*(^iAxuE`{eqGnKVc1hg&hc>n?B>^o(2f`|wbzNSgD#+oMC*ZFy^e`}0k^GhSDVj7 zXRB$tTMvW7cRV6*{*7-BpC7K5obL|rUQS;hTi;@#-2*)M_@0knzh2q|`sM_@kKeaE z%<_AW|0*T;gjeh1_f7RR>nBUc@%DIvf5Y2c94OLl;I-s__l7G|zM$9BKE^=ZgY@gu z_R{`K#_HYjb|8Mh>(|%4;})_hF)>QWgVp`+fF(iw?a_^X-h#lK<9-Icjkn1&2j&4N zF^@aT*YnAo$zjKP>xFTa2dw_g7O_N}aTvVW*H@;Rx1m9(Ja^oM5^qc`>3?r%wOOache$%2(QZ$h}R*v}tJN6CO_X`NB2Nm2Ugwlre zj2WN4JuXLKdG}`(VWV_qjMfcV7DH<0%6H*i7zxLBi6sc?E?!1{1@jZ995hzS?9HkK zoqt$D%tONVngv2}o*P=!`5j!2hQ-1Z^dU}&A=OXZi$wNW2|brzxk0YVh8PX_rU~Bn z^AxFKld20ve51x^8T?f62=iFH6hBG5M3+EVpf>X1f`)!mjve$EzjKknR4c-~yyB~E z@=4Rsa8!XDOkN-QgcwHs#DPR=+}A17Bf81-%O3OFu@WkabOv3BaE9C{Mx(q-`Bil0 zmU|&n#f_A(O)tnq!3XaPn$%6XT+m~9$1|nbhl3=72iYFmhKNj^=FceTKQ1VL1i24E ztLQ$nM^d2s<#GZvZ@s)}ggDtB(bq00*Mi(bxb9YwXRh`4&GS2!=!pFN97xFeDlQ^g z7^tEfm)yJVCLV%q_ypK2>);Etl&FJL3AWV)Y>0D?kMwghOHFK88o;Z>l>~C1CV&Zx*8?f86j#}dWmu&m^!^kee10W zIZ5lqd)N_?2vamhbTk19{j@UT9Go)sf4B?Zj$K>=g6v@OMF`f!2-hXkT$WTM)hhdQD#7P*)<5jXS9L@v@~uI+U$hekUda$B7a=+p|GJ;*0&ke( ziM9RLn(C9ZNJTkD-J=h%?(mdHK^&Voop{eu6#kP{48FVf_wnM7JNp+j-&rz$VH%jM zu@qXAwI zdUVkC6GTv6x*j#a;ZyixCVG)t<;~0jU90*NGl79d(A%$xOB0|8rD&I7?^(v_o9g~k z?i&mGo)*|Vz9@d*w|#W%4+Ihu*k+=}s(Y!u1K9Nw#{B_Hs9+Wt<`h|-^TyYH5|nr$ z6}z=jED?R!aOHerV7w%(7W+reZc!}p3I*E64}iV1{F!Fv@J zU(k1kiW7kh=oXeX#Wpdd~(~c%#O4z+OvRsqg zsCeElpy1Bc{Z&DB)CNBexRiyd)i#_AqGtRCu40B|U3z+WKoktNnygCJ2SAIQK;|Fy z2MiWZ9xyrK#l8%%EP7lSel~W@eg>C+bvW2FCsa;oSqP~|yWoaSl=se!VSO%%0yoN| zLcroZ&x3v+Afl-OvYhu8toZsswlHd^91+Ma{|X$QlT74|p_U1+-I6W8Eqv~8YnlnO zJA>51@u`WiB!MJkmj(r^3cU_|2TP@eNjaur;vwJWxx z7N6ZgItH^lln^vlFGEjUqU#@)4zknm8;h6!OrVJdcTcDtz0D-xZRW`FX&ePLPZbN* z%VP8PA>c@f1-}Fh&MSco|D-79#13D@6?LSbHJ)v%F68R~Sdvx`lKuhN>##Nm)Yi}3 zMv$di&{@_VG4Ld-o=67>B>IiXt3fn1fbfPLh~c59maB{`1R8-@F*wnzQ2%2}1|BnV z@HU?iCmk~dxC9zX&wI;CN$&_3)1fXSSC+(1fZT>3P$*!o>cox$RtqN}XHSUsp>h_$ z3l+--x>Vbfs)jAE1{|}U@WK!(DJH~ekF?ApK1Lc1dEnL852~yt?9O1dY+XgefSv@0 z;uoR(h9$Zfq6#PgNkUFuWIN8V7=l$V-MwB5ytBlU1h1)U4PpOqRWQ&pIwqNw#^g5G zP~<})pOq@%Jc_Tw&+(JSqZMW0b-B=#2ZDC3&WAJSH9o38Q&*X|YKd*vE!N(Der8r4 zlu~0E_+GHRTcjoWrq=K7=hF z3woSP#4T+Y!8{dlM7l*4MHU^=WYV~$dS%Lr0$mH&)u!9GvsG-TZ&By|>JT4OY5kF& zQ^QXtK#poqK0}B`@tbanXI7o^b z%auB0u3fb0$3t^KQ#YzQ1UrZK=d4eXj$j3)(MhV}uR9q%);Pw+(NUG8a*IhrdMEcu zOz_;27w>PP?Wi$62gHpIpwu&p+8_D|7UqamfdUrZt9zfnVJy|&qsM2;fpJrD!Fwp&>|KqOi9Wm0`MW45WA7N zT6!`}G;yFa!5x%Ja*|&f+rE(lo#~!66f4V0J&nr-3_Xs;C zCkgaW1{OKF24vM(>K~VCB%~&)1n>~d3c#rXn`P}2LxEEs;5M%DAMg<=4cUpvVj|lx zV;i>DQL<^vj>5y-ia_C|Obja$)W>m7CiCk4f#GL=!2rA=w+Tr+_OwsEQ8yxi;VB9g zkI5cBsG+W(j07ouvH@6EBpS*Zf=;Up;7n4|5mqZo*5`{Q|Uo4g&+BjNG)qvk%$cr;RsLZoo2 z>vy9fRhqzG^U+BvBSx55Y^FBgViZ_WI}!!9ua(z^dgBByx)U`xMurYi+M;*U$LDft zlw2tlYE3g|f9Ss%z*kyC;KsQ4vT(_&eF5Uo?hI>LD?!IC-#MS71zwB*y8&ZekuwGT z2~ueWn&_j_5M~KR2_?$0JVeMx*gkG|E#q1`;&0)xo5@8V$SjlRtOhH878Rn5r;mzJ z&aP}aA#LmOTrdV42e)2e4BLodbBXY!Pf-ehl1FccFr2mi)OnWzFOJVs0=%aLgLW*d zedc~xL@7c|JU4XNH!yA4*br^ns^;|*G}e#J%C+*qCcNe#C)kN$W6|<)Prgr6Z5e^V z9059?%sPKA(`2RoKvvz#G}rbChz=g&!{HL9PLsk*wlODVANwp_ak>Y(s>5BBu}Gdo zA0kB?%qgYcZ$XASW7UC^9Ofpl>vfz&c6ieao-iroRjK{a zuZTL70-Wx(%iET ziagr4Y-kN3-6IRrGVF3LP=M+W(h8em(5m`xAj!GECb*~(#M#@mVo+>Hg<+tK) zkMX7H_jvVSqT--$462}pFXj*kJqKB`;R#@6b;-XXVtdqFA!46>R2yqZ|Gpf>zzis6 zX|vupHl1k>BI#^7pT*F(GU@00zB1;4u*yZtVQ`@~=2#~^lPCX__iYw^$9H8myt#Pa zM^bJN;RYxN1(@V2o95Q`2l)1iqB^Z_w=>EOs$4-RSCHE?mi8LeR$7$Um`yH8aFl{5 z9yra1y2rYhZ5mP2!WDrUezan?A$50SSz;|$iaVG5wgtzCwo7F`NV1(dZj;4nMn}Uv zp^}z8O$jAgCR*YY;xbgSq$-;hfsYS~M--k$yyEBseY9_v62?0d_Eb7Q=oUj!%6iED zd6Rk&V#J2dZboIW9LbX$S<==#Qg_>lp{99<3dKwP&a0OtO%2~O*g9d^V?@NpxfW!$ zi{)0~^I-=GmgH-Ulm>puyB~S-toh%rzri%H99XMN?!_49bBCDHKD4gq-T+vvx3N_c zx3RG_0YQ+?%H~hhAJn$ywlP&)>`okF8pnfa!L4=Y5_9oU4Rfn*GO5x6Xka#WA$VZo ztHO4<@ur7L?bStpxai$Z$w!%rwd9j+dJ%{;NNbFKeW|sGjWd-pj$idljc?_ZaERa} z-=pplt5!S4mVb#~E|K36dGCru!o|3jB{U^dHjA81)&d+XL5cL{^qDc3w2kYzQR~9LTwSP` zUpC9$F4>J<=^goe%{G2vMqR;%IUw->q9$=JT7hoZF2s#}72Qp#blZUNyQvP(d}4V? zI$EgnSK3xa4Y;B>bz+9e=M(2Zx3D9l8A^hgLe0>mb$pk!^|gyHOe}V{IZH;L8(^*H zoT=xFCVKO~c*uR9ss$n>tuPuzjK?#ATa#hUd@PX!hec$cj?>|hDE|3AImC1%dZ^SV zz94+Z(GWb@i$7RSo#7LWDt9@FiDVQdY1=OLi3*+ROu4Bix!EucS#yyS3&YujJ<|i8 zj0Z_=tIQd2?F4m0zD{Z3S!yB41p&a*gukkOW7Yo9X4&W~Co2rlDQMJ)TtYu>olBm< z?d0KB8+6H@j3rmAlPn^}0_aRzyJyNz#`EI5Be3#0vXt|HmFGe~^}>QA3O5rco-=)+ zX9)DIPOK)Vq8HJGlyA07jOhGM+8<_1-SnlNK~-uv;0lB86N%aIo!BC*;HC!XM|)Pd z2?Ugq1(#f(>1bP+rAz+cx>mg0_(=(}Qd1(uPP;*lu#@}z9!;OGPt@13YZ z4rS99OZTBWCl!W+BFQ(n1sg~%k-9i)e$7oVRX-|Hm`q7_ec}`^pxfZy2QGd-wQ`=g z{_1vRPG725clqv>*i(2vfu22(_bs{g<>#>$E~)1*)5q36&#O8-x7T{<6$0y(;PZjs zg6H~?^*$oueWaWLf{w1;hkQR0S<#(%qP_bOf-|eBtA0vn&r;_Ch3>?%hO0%XQ>$I! zGbi=X=*c0XQfEn`rzkD}`{%QR_n4P-ko=!)fR;g&lRCN4OekN#lCego}U_^fz>D6w9b-LcdQx@%#5G-S>P{Fhy@++RnoT>G!fT0C-nXyBT8 zaqkBnLm_@D15F^vwNbM0yIIE*?;avBF_92U)@pl_(r;MWBOiI{+v}sQc5;zcpwt{@ zTmK6x&=-KJ7YC++FCXcrP&qH6ZlK#@fF}NVelWmr$$O)dZbrRLbsxG?_C%5TxKth} z$#fQdezO1|7mZS^6$F|r189MHqS8=F&zIwKbKzx&-<^!~r36mi4fNs0wAS0BulCGw ziTN|IKfEh!y_acNLK&KLf+t;PZ}xZhH+!##Cz$713QTT8KOeQ$?sMN_xs>+Lw^RiQ z;=Q!5A6p!~#_tHvb%uCnCMLV^5r=P>yrNvMo}OOsUp7yS@9G1eU-sTE+AY8FMh@eB z;e6jc{M7%@zj=Az+QGARy|%g`KA+I-@9}^zHR(w3cK-$>Dg4dD206Nmi%&oiIx_J6 zkEi9XR>m&oW~y#3)()1x`Al{T;|@7O*x{T0WozIj>&*Snnwj*Zcd$R;dYBibI9r2i zcFuU!u|D#hvkp^@-J|#!mV5hm>&F;6ua=t&q+W@qiswF=XYO4WJzu^oV;N#DuZ@d( z>hF0`3vOw2H5UkSEH`5vdhVYp8Qi|SHdQyjqB^eL0#>$j1#mKM=1&iYZ}jZ%I6V42 zqMW<+Yjj|ASKF-_V*w0x`yt9O)^JhMbZwc4LyO-Mp^>V4zv5lMnf#SdF`3Ba(eG(4R5vPm6oU((1Q7yopV{xNFJ?C?dKE+evjJ2P zJOJt${vPLoZNsu>GN|Cc8#1GyEQ*7(tc)njr&BvKFfEG1v#g9Mbuy#=>y&cu*{@!e z{sbM4cQB*?QV=cx^LZwz@tWttbz{G0GPnR-5IO+yd4?>j0O=Y3{;%e4w=XAtgwK5U z-WS&!`@NH)1<>6;ANTgb1(5>&YGxv!;57v41_=Ot2D^v6fZM=22(-ye0oyAL;-N88 zfc?iQKbg@a%px1uZEP4nnen95A{*+jQ&;@*r%Qjs>!+BBb0+J=T~FMDGAz#HvJCRc z!VJQQTCBna;kwYZcRo}HH&e$&TtbdP24MzaYcQQG=h8`cr9t^&2JPQ>?l=zSrc{Bi zS>TK9!hi67eXw1oj4=8nIio)H1LP9ZFbHRmanXEy$|&Arj)bb8zFVzc+rGqXG8tl`JIRs??}*R8E22*6&^0pP0RN%J`1fV!@mz&3z zy@|+nnMS_8teOhaYLqXhfv9DuFmDD+9v>&m!O=Y8hr+VTnsM&jtK{!ce3jp{Iuo0L2~BocMG;A{j?YYsCfvxhm;Bqd~taBTPm0pSZRe2 z(26U-Hx=*=kSs+pk3W=QRtXiM4}?Q6W4tfSGv>6t~xLvS(nVXOuZ+ zIgBhkI4L-fBqZvclq85a4mLiH9Q9}D9FN3NIiYVbCcijmf0xU%?orN?fie~#Egkt^ zecjl}>Hl2)m(Bih8%*nk2fN)s~#}9ovA4>7;NEiSo0F|djuw8b-usYk4%c0|aOoAwNWj`!0 zayJYOEMwN4lI`}A$!raSJHJEcKv*2D(E<{!lHe=DYTuF)#V*R*D$j|~3a;VEPZa*@!WYjZ#88JhXr6Q@P$PXAurNmloREJIAV9$w zfa{OH-!b@K^Y&lgf3ts3N%r3X{(T$Me<}X@t_23izievyQ}N%ocKoa2G>}vHKX3H- z6X#E&|6fROzy$IK?f*~ZKiTDfDbK?GR{o1?{-^4ntlYm;=Yc%+f8O^uHT8es=>7@# zC&~3Mz(vGA0so@B{)zG@Df2HB2jIdJ*sVV)ntuZPNs9RkfC=kwyYZVM^C!xmtKGj) za)DeNU^o6_?fWOfzpp3$(gy+YV+8^EkJZJW>VMXWzW}z`e|O;jr)vBO`tN!FFGvs& fD~^BZ(?6?#k}M?9lYebnKmf4=1_LwqudDwD(wI+W diff --git a/PSmapping.xlsx b/PSmapping.xlsx deleted file mode 100644 index 4df5af44f3cdb7440289f963da3dd128ab636601..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 12437 zcmeHtbyQqS)^FqP9ta*hxVt+9x8UyX?k)`km*DOmoZ!I;-be@%+$BihHFxfOcP4Y+ zTkrpGPOsH{PSx4JuG&(yckQht2L+7{fCaz<001(82v_iwJp=#%4FdpR0pKCui#gc4 zn%TP=sChY>x#%%^+S!s8Kts~x0U*Ke|M&V|%z-MkQTxx#=v|rzLgL*TQTY)S74Ssu zMDHp0;L<&lv?yD42xuPtsLhqt>*4#}DBDiGr@G|FuU?_j78TyIWxUsZM`?&u8dvRn zF7QkGqeO5~l}216Foq)!ALaBE_G>eoeA|jKqc!<6`5Ro0CPcO$eHV)RqLqj44;XQmnji2M@PRL#Xq zBSQxvABj0zF5Scu-Bp56*K!I(NWP_>qY^%W?d6Ha|7)KAi*52>re2k#@UEX3Is8=mF=FI)ehr8&D(4|0-9e`2A1Jem)*M$z zNx0TSO@yxYDg;U@pv(X93uui$?qG!M<~w^;3??=oMT#bE>3M)#d<( z`~3C%O@@@L7ma&Y0$pWWS-#xp28HC;bMbogDMrm#2w25L;kW`BK{`VU+Uv$YYar%D zRgbH~>sq<;4w9#`0v1zC_ECf)dE}2~GjK*+j4hUHd`E1_uJ2x{t6K6{)fwkH^OAcR znA-MUif8p;{qlXQm_4FG&W`z$bxLaZP2P2YP80k2SdP~K7easK$k}Lc3`yBKc*pur zBMB1xkah{?NITd_UIE}CJZ+i&VJ99A&Njvl4mK~l**|Rt0=%VxQT}&dRVoT{1K=ec z=_!KQBg-8JW7(C7{7C%}4SuwV7DPqH=6|(8NY}2XwIIg?3-jpUu<{rQHPZ#_lVH*aLoPP=?oXOG1QCTL zvSkrdkeW>bDb#PN_D`_UJ`OfnvG0+X{$6zm#nM(o$ z8GaNmz|7tp8U4>bi5vKf&PQ_vuSuJ}!h7T$qDnd93Jw9^a97fZ5h;=Y+9a2qf{YIJ z(@j^=hcC1}pbLhvt5fss(c$?Gu%!P}Wi+g}&ji5gk^&6?pn=~2tIWSDOO@J~J%|~- z6K_uly_3&dE84D41vr5Eij+2ZPQiO5tV9FXC|XUXqDA>?qa&$ep_++yC8o_6#lpg! zqi0Do$rr*TOoJ9>(pOg4DF>-pW`B5o!csQMO@2q-8irONuo!O{#C4g&32I(=9}M_X z(&SpB$Uk*~NoPLgih}5YXh*l}Xw8{hjZ(Lo3aCla-+TC>KNwmdAZn~1k{=`LsAajX zmZytmEX*2lWk+kxm#6;7;fTWqDj{~#n}@H*RcRB0la(1v2F1;WSHCIgLi-Lu(2unz z>sVOLxJn@EB`})o3$=k1EVyH<50CJ!>nL>o+DD_gR$7*_(LMtS&P#U6slww76l`@+ zDJi{*Qb2Q6p7yKhdmoSpaY4FO6bsgqA3vd=+4h~WO4@9Y@_1jf1}ev6{!`DvypciL zt2MHOvribXbu>E54Ltq?@0L~snKhG!Mm{hmrJ~^K5&kSlc9<+le+(d&N_IKLbHS+tO ztT~(hJvXQWC9bqN3g0KmPvP4S0lo!TWv(?R@uduX@Z z$wbxm5hsR3Ln0oK>Ya^@SzD2GF(3R<$e*3RPZcEg5Dgbv@*+unwIqtKDU@i7r5+o+ zpnJLrKq8z_3=a?+YcL(EM&5egbdM2DZ(cwB5DKK`WvROff4}VUx0>c@kkpzC*0K$J z008?Rn&x6@gsG|Z{hN?Hu+{KJFLi}hi%ZWJqJFXd|%^0P!om;)z^G$RfEZ? z=~u4wJ&P13nUGKY;grvULS5ZFp2TSLa{MTRt%KR`dFzb%r@m10g9+Q%Du|qLZerrx zePeZ^;dnm6HydTeesT4LI9@?u%zwe|#{X$C%L?QUb!Y!|vf-h|$Crz!h3od(Xt7&B za3QYh(ZiSbh}Y}RmCKgo-$<1AZ87hn5!ZQM@&+fU2S&$_kx z>*qV%S?@rDb_VI+e1t;0G9Mk~X_T&spnS(7Q(><}Wok=B&B8<5DuZY) ze_B;J!j{B0CU(rch!rVeRMH+j__{@p8MN}sm{ zEQi(yODwNLD0w+$C|n}8KkcH%T_Ez-B{9b6zR5uxaet4v)ZgI2H$5SZVs38PLf2$u z_H$mQy2ygDmfnsylO6XMi?AJ1#(f}~Hl5OUC7L$T^m)O8#8rMph>~iFgiKwrf#Sg; zE1+(|pdea+;r(IX0xxS3alsEtR!T$#$`UsL8WD2YwxpSf`>DeqXiO(N;1XG{%_3}6 zHYBFZ0)v9r@d|08Ng$9emQodAk{dR5&_+Nx81Ic^s)nfXvg}fUjaImKrhMB%*g*41 z95Qt}m3w$3;)(Two)1d_c~`bHgV6CHDMm4c5R;$%i6L_*{#_ygW7SgstbVpu3!cLa z*LKlC^C_9wY!GRm-b97X22_-59+2)QQhCZv*9Vdr$hk>vo#QBpxp-vl>L&1u-;vL6 zSEFk!#OH7uzu8gbA#U~4;Gimc4=_!0-n7Yx6uX>iz!a|`iX_%RMADEQlt7e~$`>UT zE>H;%2&sy@lSsW)fT%Zt6y9FPYK8~8wD~rP_$1n+k3~qX@9B7|f1nAHP z>Y;zF51#Fg%8SIK1o|07c~=$64;%pKibKQ8#kX=1FnQ$j6^r2)s=gxZV-u?%8CgL| zPK-`ann}x=>0&05+3?x#IJ;`%!2#1t+EQXNGAXH`!7ZC?L*5nG%G9mSbxkSMsdhVw z@4Z}0=3WzMka}rXrV`ril(g2>r!le}EfsW_kD_oVF5!rSu%(M^LI)Deif1_F1Pww5 zYe|cpKE7jkSPDj7#-ZHX5@e;+c}M9CgW!!DGQe0Upb5--K%*=uuYI|jN5IU;ixNkz z6PYpwj=RzcV_gk34sok^oEdr|A=3^n-_xg6f}#aBfg3b^zp$s4X} z45_|Xd_g76#8fWjR|2d!stkTM19+Sr9VI8-IZT`Zk9f=)qLLoZ8{cdf+c7G)GNjm}`d&+cv?otUL2rdD|mW&wa3pjqsKq zX!Mn_I}=J#Q8sEAsD*=6AEr$C3%}ogV7WdXp&TfG_x~f`s4eLX&-6p#Sn5~`!QV!MS6KWaqIG}^f*1Km_E%9_m*8Zr|z5Ozu`I3w=(B{Pa3qxs3NBDQS0{K;`;n+>rwmf+goo$E!rBvI=m; zki#b9p7JY`!9|nB$G00Dk*(s=7$x-~#za9NvZXv0iA5?bUh>JCS-+1WB=yo}{gH1TW_d>K*pGhJ z(P8Y9qi6-kGu!Y@e@!Pr8`aw+4OfK_YRZYecliRs8}p>+GBysWHs?|b_@rmgELc@8 zn7rRtp}Kgj+tYYa3y`ncx{vGLxYc~Z1|&!_&&I}_d_8Pw7hoE;&O4g<-m$nhUkR&! zX<)FmdcHl}ORu!Mq1TTF3n4xOz7M};V6mEmT>Rs9@hncZ^am$$CuCUcP=w81^pC}m z9ZO9hf(4(3=a>szoX$z7un0$0g`7LfX8n?y4Mul=G6aAUXe44niOa^?M^Ux~VDknl zHOZBdVXq?(^M~!qRe(duD)Gwh$v5So^;H&%=nFnP z?KD#F?KBQkqvhuE|J>i+>`#UhGVez2e@*nxtsig0|K`UcNuf}#@8$OyJa^xUt$<&f z$p-#UtEi!UxXPYSJ6jvMg8onQr9A<+Bao+*KZPUEti)AZQ-}%)`Gf1rY0}3>MJ}co z!pTa)VNm;82i*)pVb2aA*$8Xo#=ot}3LDz?j$QJ>F8*T9K(gMylf@_K6`0gcqh0{$ zCi+=n5KnHEcU_v+Y}p(7njZ4?)+eq!`NDpYwPOqc;XMBWclmJXXijU#13A{3Es#>nb=UpAIuHz(k1O1>G^J1*)k zosKay=E>-2>cTUC@fE5&M#7TfQ;)PHZ_P)6YGu*HeJeh~47e-l09#%OycVBbx4!M6=DO!(C=29sjV z0;x}7T0KWane-yf%XoCh?DNE}Qy*12O$x*^*BpPrdB&{2>|B>eLGFZXsYF+@(qUz@ zi|iM!&8HB^_5RDByg?r3evm$;8>d&GyN$!!5|6t_^$v9F=`-WznuU{@q2h;7))pq* z*rO8b>)g**?8L%_D;i~8JJ$4ENSiDhH27^+!6_)RZ9gNU9fbziw7zT3m$+r_)|Ap_ zo?~k~RcPBqkJ4^fKQkxX{Nx7JYg9UWWJXCvOU&q>YaG{s&D%LKN z0-JtOQ#MV?W4a+l4K-jioOw{rjZB6~jD0SkU1L&h9{)~z8^h7 zMnC6%X4Gr9UA<^7k6NtWmk{w)+-}xLGj^itO0TgVtyXM2DBj@qo3CxCg`EW?Uw~(< zxe{Nw?yWIr7*tcDbOYX#+Xu7MqJoI8%t&pRIr1pvbkf22aVzmePZh7Us@TPyDi?&& z*L9rV0#jP!h$z?OSu^BEr1z{AbxeyIYxB-BG_(TLcWcp7Q+-kQ*U^TD_W~!93Grpe zZPd&F!|Y*DdDgfSEMZs^DL|ZS4{PHuW|~4OY8)S$O*oO91Pk|!bNf;vYxc?MNX^V?&77Nb63LVJyn)_2Ph=D~a5(vB=+SPy16{AH zS@BsjjpY6qTBr@``%H+-Hy*5V**B{Kl=!%Q%QDVXglEyKICrlwa;N$8?Z6S)VVlSV zl59%JBGV_gRuyaS9PjciuWRd#N)7L|6~A;XaeVczqwi*bAIi3g20djr&GhM(Ufb?Q z1BI-bn)GH8x_W(THbXe_?m5BX7vX~<*y2BO4G5bKs%s3a=lMdCK!==i)xD>w-u}~j zDBeMfwkRijd(uH`iUN+m{2RA5L-Mxoab$Xy={Qz(s<$|KC%QqOJ~EjYZfE#z*9V@_ z>4O6P_Q1r_%%6A{9Q}raa{(y-*vwp9y==|?2u3>R!_qd&zD1k^nWjLBuqYRHth z=@rQY0y!zWd85AST*{et6vCWu@K9T2)uZ57{%p8sym6UY6R=?7*+7gxDd&_Q#g4Fm zkW@>pFm_J)KGph(Yj8x1oc?o45$(Ei$W&vi1Ez^b_Ied0xuBym* zhILseQ}b-F7W4HNuDFsS))U}2e{5(h19U2<2xh#UjX?NryTA``nv)jtRP2fr&hhi~ zLho|a49rdhMIavZw2iy-eLnYL84JnD4yH)r(m@UABv1ys^D!b=gqhpFSAgE?n0(D6 z^?uy{qaiTUx__~ z`|IPHAS!8~i{`ld28nwu`*N5vKlis*h{htN_Tcm7%FhUSE@lX&&0(<|@2l5p)lnYb zerV2XW+=zkrj=>4*#ieFfv&P&?&1jRnXaLZ?Bcd1w?FG^(uNzksVp&(CHgOKx`p?X z6u2!xE#cGj5aHS#LOr>Ziha9fFYmV!&kR3Zy=wck_mPx)RFc~NI%12@<$?k+E+N2^5}O#Tr{eMSF^=$l^4tx0#vt1 zNwl`cxGlO>UW?(QyGz(KY~I|Qu?!Ah!|Bw=1(w`)2-f?Uy^q8uSv+0)7KX*u0MR zuuK1D*o_qG+TL+mk!mZ~Ggf=|Zl?xl zed|pVZ%0zZc;i(W-{VF5Zo;RjfAY;LS#3{uuK${PL^9&J(#hCO1S)S1kvrGNj|u8I zLPzTB8zI~qUBe4K_?u`uH8q(BGr8NHm+e~d?h8IgRNSuweViA1t2QDDr8nZcn7R?I zP`^<*lArnc`8FhU4RD!&!5H9Aa-5u7AFz--{GB#nt@XKL{fXwwk2DRGtg?4`52W_Q z0n#%^bI4^u7Wslldd-?-bBm+-N3t@0!FRV$N@K$}62Gvgl58~K0O26xsOZ!gZ0pmrqSj4 zz{P5fB5fI3DLU&g5*_-B8;HVO)cVfxU0Dbj8O=V2NFEW_A_f`!hocs$@le)}p}hE* z!43LY(6Y`huQ7s!m+gH9+V1wvlt*o`c z1K@O_!b8osIT?mBV|+MrH z33AAS+&v8{^XKo`avs^CX;|>`iSfv3GWf1g45--o=u|B~$4~A@wzj!acsRdAF-!7b zz~TzbNO9&~|Kva1$~aT`)vLC(q}2y+Qf0VEU}TbG;RlA9g)4y^9`csJH|c7hgeRb{ zPf1Q|hwD1f&jm&)D~7Us+2^IOWdI#)zd0AFtt_OBosVK3I&jUw;7m1hL+n!1`)V?% z=;&sQT%(8@%eae&go38vr#Lgk28+@KMG?z0wa1OnW8Oyy`=Z3Wwp&xn=%1WY|-ye7spxNxNv$jW{iD{^}R<{X1@X_nzhHQ8l>fLQ%k zjC{Eu;~#%C<9XAzlxzZ&kMKiYT8vU#h;!16BT8l&<2YX%#Y(QZNiJ#O)+2RczyG?G zp;fTN0y{VvgfX-JU|!i zcvtcZ6pkUD+R)(BmqxLd zMkXxu>oibdR~Sq?05KoMF(_07!hIh$vS`-~r9gjn7(}&)LKC{%3xkQMeI&ANhQ^Na z>%MoXn?m7N?^2&jU?{3ZJ$z!C#56lqFmaHmC8jQEK5&p~K%Z0Ot4SD&V}Pjw89O2) zP-(;%$C#-y&q^M;3o>k)p-otnEJiP~>nOM+ze)5830t)P8+cF*cz5+_A}?U5K_YF- z^-e*dqI>ADkw+^JcA%`dpQNG@=xp$>gTZqZjbMHme*`*K1bV($9y%{Ho-|D;s!^b{ zwLU7D`NQ2FI{2{<1`}4#wsDV;yah^J<0Fh?{GIFRd5IEXwyWak)?3 zn27lnU>=Bio=yQ0=Q?%8Z8_l~w-#6t32p)P(g1`w(^Ul14iwV!3!KdG!DnMt%QnK1 zPEFEYlw6cJ=)-oUkv8fwlF97mfiWw~SJ{%miXozV>|matufJx08SZH!kPimyiWpn; zbn`_BdMq!kTuE2ekz<71!O&~)Au30jFWebx5&OO)z#hltlz*Ft?u`bXo$0c-^rJl< zW$StI9=oeF_p)P;&DIOoKxk2EwbSgx6!(aT`L!=X#N#t5bhr7wL(Yy;oMQAbNOX_< zMKZ-xKu<4N1rzcm3os_^3+6r>r9P_2*L3(Bf*uCYa{J-`f16%;!tKFS;tS6U_r zeJ~~+zudO(|S7MOqstc>vQDWFVCQoZhA6cXMQItC3zh+6u-NU25g zMfG<`SufH7Cqy;0p-5fe@iZ`8XQYKpvyUqWYD9`!hP`lO5YbE@RTRuN9Wumx9B1&9 zap^yJ@xUu2OjzOP;TQwYr9ZroYLWdqtGl#|$_y8er7oliW&HeNN|{6%X%^G$A2dSE zV?|RJ#CMAiGlw=s7DGi_ysu{0WBEDp`jB7DAtxhE(E3$lM6(pWsn09({bHlc@J%r- zmKO?lsUods=NTaq?WnGE}lvj@T*Y)g{U0XCeNqL&-}?WMXRt@s?@v z>aaiS*E%{^szk)V=Gz~Bw=F>l)cG}N4*DTsZEzE-WuwG+E z47aF1o5wvJBa}XV`pH)Nam<8m8SmY_0)aJ8#e=35gZaUa7GAg0sXg4ehzj-EZqMgq zi7x75rCL?3@wowhb@hj*81Hvtc(-yh-MLhyk7{!UUwY$5HQ6(%zp4-Lk zp!@fBSzNu(;{xNliErI2;=B9WD`t2K5JWI6XGUbxc)azPws@kmf9PDa&j1SrxSMDg znOKd|EgAShLxI;gZ8X3plh50hhh0YaeDyGcdYDAKky%4@C;c<2OOEHGa=o-O%3~eRJc}a`$6L-fJ1cVWhIcTQWjWU%U53dbV`1q#ju}Ui**FzN}+EiPxL%7>L zYF@?JfvAqF5-=`QAG^K1_bZ1|PTb!Ba5AvvwZ3)N4^qnWT)2@uy4qPYFqVZy%xHgZ+JD0g0E;0-{ z9FX2Gi9i`N&ncDa1X|73d?Rd89G>R^zhuVNZy2BkD`;{U&Al%IiAKLg8F^aq)WW?PTZ?7B_kbbMEk_#B%(U?OALC`FiyJ2qKM1Cx$D=HbBc*2Jdq zc$ng;#nzM3+4{47eg^OZK2&403&p~vlV5?8xJ&?bQl^WpGX8vfcK55-c*oCcd5}m# z4?&oCZL6qgX+jVj?ByY52YYW3ZUcndwUhFvBZs%}Ce_t+e|4#KnPdontIMQsP3qfK zD;6u`?0$4PJ+*dLfnt}u+eZtA!Gi1zKe0S^Jh~sk@pO0GdSlUAyKfZ#`8i4tR8!zZg&X$^ zDTItflC_Wm8c{%~5ntJ%Np?4UJPrsrbxQI1hBXn@cji|h@3wwf6PkpM4k;@9=LiGA z56>t;e_xmZ!A+1D1714upa1~$e_c9_93B6=a)K@S&nr7g!2$GA^#yf90=-7avj!I~ ztie)5@4Zw9knyoFi^4zbX~f8u2Ks7dZPZaU|3fdxmlNyj+S|8 zsjaSdr~l;SAO<^(l>M7YuTE~=2;t>Z&iZybb8cOqa$~e55$!r6Ei03g zK<)`29MFPxXWY%Fe&K`m3PE~)%y8JON>|Pg$y*hrt(lxNsJ=T+9e{-4`DuFRjs$2~2mckEjfk=y)5vcZV z!fvEF;Yv8xiKeFE{x;VTaYTx4f?}u2IwXkmWtZ`!3o8zb_&&W*W0ix}gj*L@mI#^Z6)q!>M?8P($75L5UvNG=!^8wsdjyu*j5$#8Xg?~SSPh05MY!u7b>^wH1RdQIjy)fJ~ovJsy%iz zW@^+D+C70q@sFb<2*|hK+K#_3t^TLL{`vbih1N=P{|@l)MW6p9eEBT}TgAVWfc{SS z?}c#xN;nQqUH;$8;(o{Zy*lbIq-XGn@o%+KzZ3sn@bVY271AHX|55((JHqeT=D!da z(f{Y~|6jT1-vNJ5-u?y1jPVEH@2TA10e(L){|n$3-aqE_KgZ|4lm33R_!lWE*mwMs zm;Z6l_&dtKdz61+004S#0D%ASFMlWh-5vM~U<2$O{C_$Jzk~k!QvDYs*z2JCmofdh Ya4X5dfOQZ6Kn8!Pz-p}h_T}mS0H#OTDF6Tf From 94b7c0e8d5854d0dfaa66820685e7c058e14fc58 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Tue, 22 Nov 2022 16:42:58 +0100 Subject: [PATCH 06/23] adding a script to display LGEN values --- DisplayLGENs.m | 90 ++++++++++++++++++++++++++++++++++++ operations/LGENvisualCheck.m | 64 +++++++++++++++++-------- 2 files changed, 136 insertions(+), 18 deletions(-) create mode 100644 DisplayLGENs.m diff --git a/DisplayLGENs.m b/DisplayLGENs.m new file mode 100644 index 0000000..93375bd --- /dev/null +++ b/DisplayLGENs.m @@ -0,0 +1,90 @@ +% {}~ + +%% description +% this is a script which displays LGEN currents; + +%% include libraries +% - include Matlab libraries +pathToLibrary="./"; +addpath(genpath(pathToLibrary)); +pathToLibrary="../MachineRefs"; +addpath(genpath(pathToLibrary)); + +%% settings + +% ------------------------------------------------------------------------- +% USER's input data +machine=["ISO1" "ISO2" "ISO3" "ISO4" "ISO1" "ISO2" "ISO3" "ISO4"]'; +beamPart=["PROTON" "PROTON" "PROTON" "PROTON" "CARBON" "CARBON" "CARBON" "CARBON" ]; +% machine="ISO2"; +% beamPart=["PROTON" "CARBON" ]; +config="TM"; % select configuration: TM, RFKO +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% check of user input data +nSets=max([length(beamPart) length(machine) length(config)]); +beamPart=MyCheck(beamPart,nSets,"beamPart"); +machine=MyCheck(machine,nSets,"machine"); +config=MyCheck(config,nSets,"config"); + +%% clear variables +[cyCodes,ranges,Eks,Brhos,currents,fields,kicks,psNames,FileNameCurrents,magNames]=... + deal(missing(),missing(),missing(),missing(),missing(),missing(),missing(),missing(),missing(),missing()); +clear PSmapping; +myLeg=strings(nSets,1); + +%% parse DBs + +% - get PS mapping +FullFileName=ReturnDefFile("PSmapping"); PSmapping=readtable(FullFileName); + +% - get values +for iSet=1:nSets + clear tmpCyCodes tmpRanges tmpEks tmpBrhos tmpCurrents tmpFields tmpKicks tmpPsNames tmpFileNameCurrents tmpMagNames ; + myConfig=sprintf("%s,%s,%s",machine(iSet),beamPart(iSet),config(iSet)); + [tmpCyCodes,tmpRanges,tmpEks,tmpBrhos,tmpCurrents,tmpFields,tmpKicks,tmpPsNames,tmpFileNameCurrents]=AcquireLGENValues(myConfig); + tmpPsNames=string(tmpPsNames); + tmpCyCodes=upper(string(tmpCyCodes)); + tmpMagNames=MagNames2LGENnames(tmpPsNames,true,PSmapping); + myLeg(iSet)=myConfig; % comment me, if single data set + % - store data + cyCodes=ExpandMat(cyCodes,tmpCyCodes); + ranges=ExpandMat(ranges,tmpRanges); + Eks=ExpandMat(Eks,tmpEks); + Brhos=ExpandMat(Brhos,tmpBrhos); + currents=ExpandMat(currents,tmpCurrents); + fields=ExpandMat(fields,tmpFields); + kicks=ExpandMat(kicks,tmpKicks); + psNames=ExpandMat(psNames,tmpPsNames); + FileNameCurrents=ExpandMat(FileNameCurrents,tmpFileNameCurrents); + magNames=ExpandMat(magNames,tmpMagNames); +end + +% - normalise currents to Brho +clear normCurrents; normCurrents=NaN(size(currents)); +for iSet=1:nSets + normCurrents(:,:,iSet)=currents(:,:,iSet)./Brhos(:,iSet); +end + +%% visual checks +LGENvisualCheck(psNames,Eks ,"Ek [MeV/u]",currents,"I [A]",magNames,"",myLeg); +LGENvisualCheck(psNames,Eks ,"Ek [MeV/u]",normCurrents,"I/B\rho [A/Tm]",magNames,"",myLeg); +LGENvisualCheck(psNames,ranges,"range [mm]",currents,"I [A]",magNames,"",myLeg); +LGENvisualCheck(psNames,ranges,"range [mm]",normCurrents,"I/B\rho [A/Tm]",magNames,"",myLeg); +LGENvisualCheck(psNames,Brhos ,"B\rho [Tm]",currents,"I [A]",magNames,"",myLeg); +LGENvisualCheck(psNames,Brhos ,"B\rho [Tm]",normCurrents,"I/B\rho [A/Tm]",magNames,"",myLeg); + +%% local functions +function OutVar=MyCheck(InVar,nSets,myName) + OutVar=InVar; + nIn=length(OutVar); + if (nIn Date: Tue, 22 Nov 2022 16:44:30 +0100 Subject: [PATCH 07/23] removing unused declaration in CreateQAfiles.m --- CreateQAfiles.m | 3 --- 1 file changed, 3 deletions(-) diff --git a/CreateQAfiles.m b/CreateQAfiles.m index cb1d173..40e0c14 100644 --- a/CreateQAfiles.m +++ b/CreateQAfiles.m @@ -19,9 +19,6 @@ % ------------------------------------------------------------------------- % USER's input data -kPath="P:\Accelerating-System\Accelerator-data"; -% kPath="K:"; - machine="ISO3"; beamPart="CARBON"; config="TM"; % select configuration: TM, RFKO From c0e2323c203e83ad0d18d3b70e5e84c5692c4682 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 23 Nov 2022 11:53:16 +0100 Subject: [PATCH 08/23] Adding a script which shows basic beam properties changing with range --- DisplayLGENs.m | 2 +- DisplayRefValues.m | 140 +++++++++++++++++++++++++ educational/PlotBasicQuantitiesPandC.m | 110 ------------------- 3 files changed, 141 insertions(+), 111 deletions(-) create mode 100644 DisplayRefValues.m delete mode 100644 educational/PlotBasicQuantitiesPandC.m diff --git a/DisplayLGENs.m b/DisplayLGENs.m index 93375bd..12041d8 100644 --- a/DisplayLGENs.m +++ b/DisplayLGENs.m @@ -14,7 +14,7 @@ % ------------------------------------------------------------------------- % USER's input data -machine=["ISO1" "ISO2" "ISO3" "ISO4" "ISO1" "ISO2" "ISO3" "ISO4"]'; +machine=["ISO1" "ISO2" "ISO3" "ISO4" "ISO1" "ISO2" "ISO3" "ISO4"]; beamPart=["PROTON" "PROTON" "PROTON" "PROTON" "CARBON" "CARBON" "CARBON" "CARBON" ]; % machine="ISO2"; % beamPart=["PROTON" "CARBON" ]; diff --git a/DisplayRefValues.m b/DisplayRefValues.m new file mode 100644 index 0000000..a19cb94 --- /dev/null +++ b/DisplayRefValues.m @@ -0,0 +1,140 @@ +% {}~ + +%% description +% this is a script which displays basic beam properties which depend on beam energy; + +%% include libraries +% - include Matlab libraries +pathToLibrary="./"; +addpath(genpath(pathToLibrary)); +pathToLibrary="../MachineRefs"; +addpath(genpath(pathToLibrary)); + +%% settings + +% ------------------------------------------------------------------------- +% USER's input data +% what to load +machine="ISO1"; +beamPart=["PROTON" "CARBON"]; +% what to show +showX="range"; +showXlabel="R [mm]"; +showY=[ "beta" "gamma" "betagamma" "pc" "Brho" "Ek" ]; +showYlabel=["\beta_{rel} []" "\gamma_{rel} []" "\beta_{rel}\gamma_{rel} []" "pc [MeV/c/u]" "B\rho [Tm]" "E_k [MeV/u]"]; +nCoupled=3; +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% check of user input data +nSets=max([length(beamPart) length(machine)]); +beamPart=MyCheck(beamPart,nSets,"beamPart"); +machine=MyCheck(machine,nSets,"machine"); +mSets=max([length(showX) length(showXlabel) length(showY) length(showYlabel)]); +showX=MyCheck(showX,mSets,"showX"); +showXlabel=MyCheck(showXlabel,mSets,"showXlabel"); +showY=MyCheck(showY,mSets,"showY"); +showYlabel=MyCheck(showYlabel,mSets,"showYlabel"); + +%% clear variables +xVals=missing(); +yVals=missing(); + +%% parse DBs + +% - get values +myLeg=strings(nSets,1); +for iSet=1:nSets + clear EnData tmpDataX tmpDataY; + myConfig=sprintf("%s,%s",machine(iSet),beamPart(iSet)); + FullFileName=ReturnDefFile("BRHO",myConfig); EnData=readtable(FullFileName); + myLeg(iSet)=myConfig; + % - store data + [tmpDataX,tmpDataY]=ExtractData(EnData,showX,showY); + xVals=ExpandMat(xVals,tmpDataX); + yVals=ExpandMat(yVals,tmpDataY); +end +xVals=permute(xVals,[1 3 2]); % what to show is the outermost dimension +yVals=permute(yVals,[1 3 2]); + +%% show stuff +ShowStuff(xVals,yVals,showXlabel,showYlabel,myLeg,nCoupled); + +%% local functions +function ShowStuff(xVals,yVals,showXlabel,showYlabel,myLeg,nCoupled,myTitle) + if (~exist("nCoupled","var")), nCoupled=1; end + if (~exist("myTitle","var")), myTitle=missing(); end + nShows=size(xVals,3); + nSets=size(xVals,2); + + figure(); + cm=colormap(parula(nSets+1)); + % - set grid of plots + nPlots=nShows; + if (~ismissing(myLeg)), nPlots=nPlots+1; end + [nRows,nCols]=GetNrowsNcols(nPlots,nCoupled); + % - set markers + if ( nSets==2 ) + markers=[ "o" "*" ]; + else + markers=strings(nSets,1); + markers(:)="."; + end + % - actually plot + for iPlot=1:nShows + subplot(nRows,nCols,iPlot); + lFirst=true; + for iSet=1:nSets + if (lFirst), lFirst=false; else, hold on; end + plot(xVals(:,iSet,iPlot),yVals(:,iSet,iPlot),"-","Marker",markers(iSet),"Color",cm(iSet,:)); + end + xlabel(showXlabel(iPlot)); ylabel(showYlabel(iPlot)); grid on; + end + % - legend plot + if (~ismissing(myLeg)) + subplot(nRows,nCols,nPlots); + lFirst=true; + for iSet=1:nSets + if (lFirst), lFirst=false; else, hold on; end + plot(NaN(),NaN(),"-","Marker",markers(iSet),"Color",cm(iSet,:)); + end + legend(myLeg); + end + % - general stuff + if (~ismissing(myTitle)), sgtitle(myTitle); end +end + +function [DataX,DataY]=ExtractData(EnData,showX,showY) + names=fieldnames(EnData); + DataX=missing(); DataY=missing(); + for ii=1:length(showX) + isField=strcmpi(showX(ii),names); + if any(isField) + DataX=ExpandMat(DataX,EnData.(names{isField})); + else + error("no field %s in table!",showX(ii)); + end + end + for ii=1:length(showY) + isField=strcmpi(showY(ii),names); + if any(isField) + DataY=ExpandMat(DataY,EnData.(names{isField})); + else + error("no field %s in table!",showY(ii)); + end + end + +end + +function OutVar=MyCheck(InVar,nSets,myName) + OutVar=InVar; + nIn=length(OutVar); + if (nIn Date: Tue, 31 Jan 2023 14:44:09 +0100 Subject: [PATCH 09/23] Adding a title to the window of LGEN plots --- DisplayLGENs.m | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/DisplayLGENs.m b/DisplayLGENs.m index 12041d8..5a064b3 100644 --- a/DisplayLGENs.m +++ b/DisplayLGENs.m @@ -14,11 +14,13 @@ % ------------------------------------------------------------------------- % USER's input data -machine=["ISO1" "ISO2" "ISO3" "ISO4" "ISO1" "ISO2" "ISO3" "ISO4"]; -beamPart=["PROTON" "PROTON" "PROTON" "PROTON" "CARBON" "CARBON" "CARBON" "CARBON" ]; +machine=["ISO1" "ISO2" "ISO3" "ISO4" ]; % "ISO1" "ISO2" "ISO3" "ISO4"]; +% beamPart=["PROTON" "PROTON" "PROTON" "PROTON" ]; %"CARBON" "CARBON" "CARBON" "CARBON" ]; +beamPart=["CARBON" "CARBON" "CARBON" "CARBON" ]; % machine="ISO2"; % beamPart=["PROTON" "CARBON" ]; config="TM"; % select configuration: TM, RFKO +myTitle="LGEN values as downloaded on 2023-01-30"; % ------------------------------------------------------------------------- % ------------------------------------------------------------------------- @@ -68,12 +70,12 @@ end %% visual checks -LGENvisualCheck(psNames,Eks ,"Ek [MeV/u]",currents,"I [A]",magNames,"",myLeg); -LGENvisualCheck(psNames,Eks ,"Ek [MeV/u]",normCurrents,"I/B\rho [A/Tm]",magNames,"",myLeg); -LGENvisualCheck(psNames,ranges,"range [mm]",currents,"I [A]",magNames,"",myLeg); -LGENvisualCheck(psNames,ranges,"range [mm]",normCurrents,"I/B\rho [A/Tm]",magNames,"",myLeg); -LGENvisualCheck(psNames,Brhos ,"B\rho [Tm]",currents,"I [A]",magNames,"",myLeg); -LGENvisualCheck(psNames,Brhos ,"B\rho [Tm]",normCurrents,"I/B\rho [A/Tm]",magNames,"",myLeg); +LGENvisualCheck(psNames,Eks ,"Ek [MeV/u]",currents,"I [A]",magNames,myTitle,myLeg); +LGENvisualCheck(psNames,Eks ,"Ek [MeV/u]",normCurrents,"I/B\rho [A/Tm]",magNames,myTitle,myLeg); +LGENvisualCheck(psNames,ranges,"range [mm]",currents,"I [A]",magNames,myTitle,myLeg); +LGENvisualCheck(psNames,ranges,"range [mm]",normCurrents,"I/B\rho [A/Tm]",magNames,myTitle,myLeg); +LGENvisualCheck(psNames,Brhos ,"B\rho [Tm]",currents,"I [A]",magNames,myTitle,myLeg); +LGENvisualCheck(psNames,Brhos ,"B\rho [Tm]",normCurrents,"I/B\rho [A/Tm]",magNames,myTitle,myLeg); %% local functions function OutVar=MyCheck(InVar,nSets,myName) From eb54ebe93c623beed36f75c33f532a2c670d6ddd Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Tue, 31 Jan 2023 14:45:26 +0100 Subject: [PATCH 10/23] bug fix in comparing PSnames to write with all PS names + repeating a single point in the scan N times --- CreateQAfiles.m | 32 +++++++++++++++++--------------- 1 file changed, 17 insertions(+), 15 deletions(-) diff --git a/CreateQAfiles.m b/CreateQAfiles.m index 40e0c14..008edb5 100644 --- a/CreateQAfiles.m +++ b/CreateQAfiles.m @@ -20,7 +20,7 @@ % ------------------------------------------------------------------------- % USER's input data machine="ISO3"; -beamPart="CARBON"; +beamPart="PROTON"; config="TM"; % select configuration: TM, RFKO % ------------------------------------------------------------------------- @@ -46,20 +46,22 @@ % ------------------------------------------------------------------------- % USER's input data % wrMagnetNames=[ "H2-012A-QUE" "H2-016A-QUE" "H2-022A-QUE" "HE-025A-QUE" ]; -wrMagnetNames=[ "H2-012A-QUE" "H2-016A-QUE" "H2-022A-QUE" "HE-025A-QUE" ]; +% wrMagnetNames=[ "H2-012A-QUE" "H2-016A-QUE" "H2-022A-QUE" "HE-025A-QUE" ]; % wrMagnetNames=[ "HE-018A-QUE" "HE-020A-QUE" "HE-023A-QUE" "HE-025A-QUE" ]; -wrRange=[30 30 30 30]; % [mm] +wrMagnetNames=[ "HE-H07A-CEB" "HE-V07A-CEB" "HE-H27A-CEB" "HE-V27A-CEB" ]; +wrRange=[320 320 320 320]; % [mm] wrScan=["scanTM" "scanTM" "scanTM" "scanTM"]; -wrDImin=[20 20 20 20 ]; % [A] -wrDImax=[20 20 20 20 ]; % [A] -wrDIdel=[1 1 1 1 ]; % [A] -wrNtimes=[ 1 1 1 1]; -wrIbef=[5 5 5 5 ]; % [A] -wrNIbef=[3 3 3 3]; -wrIaft=[350 350 350 350 ]; % [A] +wrDImin=[210 140 140 140 ]; % min current [A] +wrDImax=[ 70 140 140 140 ]; % max current [A] +wrDIdel=[ 70 70 70 70 ]; % delta current [A] +wrNtimes=[ 1 1 1 1 ]; % repeat scan N times +wrNpoints=[ 5 5 5 5 ]; % repeat each point N times +wrIbef=[ -150 -150 -150 -150 ]; % [A] +wrNIbef=[ 2 2 2 2 ]; +wrIaft=[150 150 150 150 ]; % [A] wrNIaft=[2 2 2 2]; -wrImin=[ 5 5 5 5 ]; -wrImax=[ 350 350 350 350]; +wrImin=[ -150 -150 -150 -150 ]; +wrImax=[ 150 150 150 150 ]; oFileName="test.xlsx"; % ------------------------------------------------------------------------- @@ -70,8 +72,8 @@ % echo TM values rTM=find(rangesTM==wrRange(ii)); if ( isempty(rTM) ), error("Range %d mm not available in TM table!",wrRange(ii)); end - pTM=find(strcmpi(psNamesTM,string(PSmapping.LGEN(jj)))); - if ( isempty(pTM) ), error("LGEN name %s not found in TM table!",PSmapping.LGEN(jj)); end + pTM=find(strcmpi(psNamesTM,wrPSnames(ii))); + if ( isempty(pTM) ), error("LGEN name %s not found in TM table!",wrPSnames(ii)); end warning("...TM value of %s (aka %s) for %s at %d mm: %f A;",wrPSnames(ii),wrMagnetNames(ii),beamPart,wrRange(ii),currentsTM(rTM,pTM)); % array characteristics switch upper(wrScan(ii)) @@ -86,7 +88,7 @@ Imax=currentsTM(rTM,pTM)+wrDImax(ii); Idel=wrDIdel(ii); end - tmpScan=(Imin:Idel:Imax)'; + tmpScan=repelem(Imin:Idel:Imax,wrNpoints(ii))'; tmpScan=CorrectRange(tmpScan,wrImin(ii),wrImax(ii)); tmpScan=RepeatScan(tmpScan,wrNtimes(ii)); tmpScan=DecorateScan(tmpScan,wrIbef(ii),wrIaft(ii),wrNIbef(ii),wrNIaft(ii)); From e5baf04927de16e8d37d82d24809dad65807a46b Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 23 Feb 2023 11:10:48 +0100 Subject: [PATCH 11/23] Minor updates to operations/LGENvisualCheck.m: * better palette; * possibility to filter on (HEBT) families; * filtering on magnet names (no longer PS/LGEN names); * showing plots in order of magnet name (no longer PS/LGEN names); --- operations/LGENvisualCheck.m | 30 ++++++++++++++++++++++-------- 1 file changed, 22 insertions(+), 8 deletions(-) diff --git a/operations/LGENvisualCheck.m b/operations/LGENvisualCheck.m index ef66393..93359a7 100644 --- a/operations/LGENvisualCheck.m +++ b/operations/LGENvisualCheck.m @@ -1,12 +1,26 @@ -function LGENvisualCheck(psNames,xVals,xLab,yVals,yLab,magNames,myTitle,myLegend,LGENs) +function LGENvisualCheck(psNames,xVals,xLab,yVals,yLab,magNames,myTitle,myLegend,magNames2show) if (~exist("myTitle","var") || strlength(myTitle)==0 ), myTitle=missing(); end if (~exist("myLegend","var")), myLegend=missing(); end - if (~exist("LGENs","var")), LGENs=unique(psNames); end % default: do not filter - nMagnets=size(LGENs,1); - nSets=size(psNames,2); + if (exist("magNames2show","var")) + switch upper(magNames2show) + case {"QUE","QUES","QUADRUPOLES"} + magNames2show=unique(magNames(contains(magNames,"QUE"))); + if (~ismissing(myTitle)), myTitle=strcat(myTitle," - only QUEs"); end + case {"CEB","CEBS","CORRECTORS"} + magNames2show=unique(magNames(contains(magNames,"CEB"))); + if (~ismissing(myTitle)), myTitle=strcat(myTitle," - only CEBs"); end + case {"DIP","DIPS","DIPOLES"} + magNames2show=unique(magNames(contains(magNames,"SWH")|contains(magNames,"MBS")|contains(magNames,"SW2"))); + if (~ismissing(myTitle)), myTitle=strcat(myTitle," - only MAIN DIPOLEs"); end + end + else + magNames2show=unique(magNames(~ismissing(magNames))); % default: do not filter + end + nMagnets=size(magNames2show,1); + nSets=size(magNames,2); figure(); - cm=colormap(parula(nSets+1)); + cm=colormap(turbo(nSets)); % - set grid of plots nPlots=nMagnets; if (~ismissing(myLegend)), nPlots=nPlots+1; end @@ -23,9 +37,9 @@ function LGENvisualCheck(psNames,xVals,xLab,yVals,yLab,magNames,myTitle,myLegend subplot(nRows,nCols,iMagnet); lFirst=true; for iSet=1:nSets - jj=find(strcmp(psNames(:,iSet),LGENs(iMagnet))); + jj=find(strcmp(magNames(:,iSet),magNames2show(iMagnet))); if ( isempty(jj) ) - warning("...unable to find LGEN %s in list of PS names!",LGENs(iMagnet)); + warning("...unable to find LGEN %s in list of PS names!",magNames2show(iMagnet)); continue end if (lFirst), lFirst=false; else, hold on; end @@ -33,7 +47,7 @@ function LGENvisualCheck(psNames,xVals,xLab,yVals,yLab,magNames,myTitle,myLegend plot(xVals(IDs,iSet),yVals(IDs,jj,iSet),"-","Marker",markers(iSet),"Color",cm(iSet,:)); end xlabel(xLab); ylabel(yLab); grid on; - title(sprintf("%s (aka %s)",magNames(jj,iSet),LGENs(iMagnet))); + title(sprintf("%s (aka %s)",magNames2show(iMagnet),psNames(jj,iSet))); end % - legend plot if (~ismissing(myLegend)) From 4a29ac8928e154bef9597309e45e86a985b4fb40 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 11 Oct 2023 11:59:47 +0200 Subject: [PATCH 12/23] adding Helium to particle DB --- educational/particleData.m | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/educational/particleData.m b/educational/particleData.m index be2be55..aabbd09 100644 --- a/educational/particleData.m +++ b/educational/particleData.m @@ -7,6 +7,11 @@ Zp=1; mp=938.27208816; % [MeV/c2] +% helium data +AHe=4; +ZHe=2; +mHe=931.2386834*AHe; % [MeV/c2] + % carbon data AC=12; ZC=6; From c599b006c03bb117d968932772a9fd0b697e67d3 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 11 Oct 2023 12:00:06 +0200 Subject: [PATCH 13/23] cleared ShowMCS.m of particle-dependent variable names --- educational/ShowMCS.m | 55 ++++++++++++++++++++++++------------------- 1 file changed, 31 insertions(+), 24 deletions(-) diff --git a/educational/ShowMCS.m b/educational/ShowMCS.m index 401c88d..eef4284 100644 --- a/educational/ShowMCS.m +++ b/educational/ShowMCS.m @@ -13,51 +13,58 @@ addpath(genpath(pathToLibrary)); %% settings -% - proton energies -Ek_p=1:1:230; % [MeV] -% Ek_p=228.57; -% - carbon energies -Ek_C=1:1:400; % [MeV/A] -% Ek_C=398.84; % [MeV/A] + +% - particle +myPart="PROTON"; % available: "PROTON", "CARBON", "HELIUM" + +% - kinetic energies +Ek=1:1:250; % [MeV] % proton energies +% Ek=1:1:400; % [MeV/A] % carbon energies +% Ek_p=228.57; % single energy + % - range traversed -mmEquiv=1.0; % [mm] +mmEquiv=100.0; % [mm] % mmEquiv=0.1:0.1:10.1; % [mm] + % - water material parameters X0_H2O=36.08; % [cm] -%% compute MCS scattering tables +%% Load particle data -% - load particle data run(".\particleData.m"); +switch upper(myPart) + case {"P","PROT","PROTON"} + myM=mp; myEk=Ek; myZ=Zp; unitEk="MeV"; + case {"HE","HELIUM","ALPHA","ALFA"} + myM=mHe; myEk=Ek*AHe; myZ=ZHe; unitEk="MeV/u"; + case {"C","CARB","CARBON"} + myM=mC; myEk=Ek*AC; myZ=ZC; unitEk="MeV/u"; + otherwise + error("Unknown particle %s!",myPart); +end +%% compute MCS scattering tables % - relativistic quantities -[beta_p,gamma_p,betagamma_p]=ComputeRelativisticQuantities(Ek_p,mp); % [], [], [] -[beta_C,gamma_C,betagamma_C]=ComputeRelativisticQuantities(Ek_C*AC,mC); % [], [], [] - +[myBeta,myGamma,myBetaGamma]=ComputeRelativisticQuantities(myEk,myM); % [], [], [] % - MCS theta0 (Moliere's theory) -theta0_p=ComputeTheta0(beta_p,betagamma_p*mp,Zp,mmEquiv,X0_H2O); -theta0_C=ComputeTheta0(beta_C,betagamma_C*mC,ZC,mmEquiv,X0_H2O); +theta0=ComputeTheta0(myBeta,myBetaGamma*myM,myZ,mmEquiv,X0_H2O); % - show theta0 -if (length(mmEquiv)==1 && (length(Ek_p)>1 || length(Ek_C)>1)) +if (length(mmEquiv)==1 && length(Ek)>1) % as a function of energy - ShowMe(theta0_p*1000,Ek_p,"\theta [mrad]","E_k [MeV]",sprintf("MCS for PROTONs after traversing %g mm of water equivalent",mmEquiv)); set(gca, 'YScale', 'log'); - ShowMe(theta0_C*1000,Ek_C,"\theta [mrad]","E_k [MeV/u]",sprintf("MCS for CARBON ions after traversing %g mm of water equivalent",mmEquiv)); set(gca, 'YScale', 'log'); + ShowMe(theta0*1000,Ek,"\theta [mrad]",sprintf("E_k [%s]",unitEk),sprintf("MCS for %s after traversing %g mm of water equivalent",myPart,mmEquiv)); set(gca, 'YScale', 'log'); else % as a function of material thickness - ShowMe(theta0_p*1000,mmEquiv,"\theta [mrad]","H_2O Range [mm]",sprintf("MCS for PROTONs of %g MeV",Ek_p)); set(gca, 'YScale', 'log'); - ShowMe(theta0_C*1000,mmEquiv,"\theta [mrad]","H_2O Range [mm]",sprintf("MCS for CARBON ions of %g MeV/u",Ek_C)); set(gca, 'YScale', 'log'); + ShowMe(theta0*1000,mmEquiv,"\theta [mrad]","H_2O Range [mm]",sprintf("MCS for %s of %g %s",myPart,myEk,unitEk)); set(gca, 'YScale', 'log'); end %% show effect of scattering L=0.15; % distance travelled by scattered bea, [m] -if (length(mmEquiv)==1 && (length(Ek_p)>1 || length(Ek_C)>1)) +if (length(mmEquiv)==1 && length(Ek)>1) % as a function of energy - ShowMe(theta0_p*L*1000,Ek_p,"\theta\timesL [mm]","E_k [MeV]",sprintf("MCS for PROTONs after traversing %g mm of water equivalent",mmEquiv)); set(gca, 'YScale', 'log'); - ShowMe(theta0_C*L*1000,Ek_C,"\theta\timesL [mm]","E_k [MeV/u]",sprintf("MCS for CARBON ions after traversing %g mm of water equivalent",mmEquiv)); set(gca, 'YScale', 'log'); + ShowMe(theta0*L*1000,Ek,"\theta\timesL [mm]",sprintf("E_k [%s]",unitEk),sprintf("MCS for %s after traversing %g mm of water equivalent",myPart,mmEquiv)); set(gca, 'YScale', 'log'); else % as a function of material thickness - ShowMe(theta0_p*1000,mmEquiv,"\theta\timesL [mm]","H_2O Range [mm]",sprintf("MCS for PROTONs of %g MeV",Ek_p)); set(gca, 'YScale', 'log'); - ShowMe(theta0_C*1000,mmEquiv,"\theta\timesL [mm]","H_2O Range [mm]",sprintf("MCS for CARBON ions of %g MeV/u",Ek_C)); set(gca, 'YScale', 'log'); + ShowMe(theta0*1000,mmEquiv,"\theta\timesL [mm]","H_2O Range [mm]",sprintf("MCS for %s of %g %s",myPart,Ek_p,unitEk)); set(gca, 'YScale', 'log'); end From d90f0c37f7d3bd3bf21eba5dd69f1e26e42944e9 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 2 Nov 2023 16:04:04 +0100 Subject: [PATCH 14/23] showing MCS and Bethe-Bloch/Landau-Vavilov also for He --- educational/ShowBetheBlochLandauVavilov.m | 151 +++++++++------------- educational/ShowMCS.m | 23 +--- educational/particleData.m | 7 + educational/setParticle.m | 17 +++ 4 files changed, 94 insertions(+), 104 deletions(-) create mode 100644 educational/setParticle.m diff --git a/educational/ShowBetheBlochLandauVavilov.m b/educational/ShowBetheBlochLandauVavilov.m index cabc707..ce0c8b8 100644 --- a/educational/ShowBetheBlochLandauVavilov.m +++ b/educational/ShowBetheBlochLandauVavilov.m @@ -1,30 +1,42 @@ % {}~ %% description -% this is a script which shows the Bethe-Bloch for the accelerated particles -% Refs: -% - main formula: PDG, 2018, Chap. 33, pag. 447 (eq. 33.5) -% - material parameters: -% . https://pdg.lbl.gov/2022/AtomicNuclearProperties/index.html -% - density effect: Sternheimer, Berger and Seltzer, ATOMIC DATA AND NUCLEAR DATA TABLES 30,26 l-27 1 (1984) +% this is a script which shows: +% * the Bethe-Bloch for the accelerated particles. Refs: +% - main formula: PDG, 2018, Chap. 33, pag. 447 (eq. 33.5) +% - material parameters: +% . https://pdg.lbl.gov/2022/AtomicNuclearProperties/index.html +% - density effect: Sternheimer, Berger and Seltzer, ATOMIC DATA AND NUCLEAR DATA TABLES 30,26 l-27 1 (1984) +% * the Landau-Vavilov for the accelerated particles. Refs: +% - main formula: PDG, 2018, Chap. 33, pag. 450 (eq. 33.11) +% - material parameters: +% . https://pdg.lbl.gov/2022/AtomicNuclearProperties/index.html %% include libraries % - include Matlab libraries pathToLibrary="..\"; addpath(genpath(pathToLibrary)); +%% clean/close +clear all; +close all; + %% settings -% - proton energies -Ek_p=1:1:230; % [MeV] -EkSel_p=228.57; % [MeV] -if (~exist("EkSel_p","var")),EkSel_p=max(Ek_p); end -% - carbon energies -Ek_C=1:1:400; % [MeV/A] -EkSel_C=398.84; % [MeV/A] -if (~exist("EkSel_C","var")),EkSel_C=max(Ek_C); end + +% - particle +myPart="PROTON"; % available: "PROTON", "CARBON", "HELIUM" + +% NOTA BENE: please choose either an array of values of Ek and a single value +% for the range traversed or the other way around; +% - kinetic energies +% Ek=1:1:1000; % [MeV] % proton energies +% Ek=1:1:400; % [MeV/A] % carbon energies +Ek=228.57; % single energy + % - range traversed -% mmEquiv=1; % [mm] -mmEquiv=0.1:0.01:4.1; % [mm] +% mmEquiv=100.0; % [mm] +mmEquiv=0.1:0.1:10.1; % [mm] + % - water material parameters ZoA_H2O=0.555087; % [] I_H2O=79.7; % [eV] @@ -37,110 +49,73 @@ densEff_m_H2O=3.4773; % [] densEff_d0_H2O=0.0; % [] -%% Bethe-Bloch +%% complement user's input +if (length(Ek)>1 && length(mmEquiv)>1), error("Please choose to scan either energy or range!"); end +if (length(Ek)==1), Ek=1:1:Ek; end -% - load particle data -run(".\particleData.m"); +%% Load particle data +% returns: myM [MeV/c2], myEk [MeV], myZ [], unitEk ("MeV" for protons, "MeV/u" for others); +run(".\setParticle.m"); + +%% Bethe-Bloch % - relativistic quantities -[beta_p,gamma_p,betagamma_p]=ComputeRelativisticQuantities(Ek_p,mp); % [], [], [] -[beta_C,gamma_C,betagamma_C]=ComputeRelativisticQuantities(Ek_C*AC,mC); % [], [], [] +[myBeta,myGamma,myBetaGamma]=ComputeRelativisticQuantities(myEk,myM); % [], [], [] % - Wmax -Wmax_p=ComputeWmax(betagamma_p,gamma_p,mp); % [MeV] -Wmax_C=ComputeWmax(betagamma_C,gamma_C,mC); % [MeV] +Wmax=ComputeWmax(myBetaGamma,myGamma,myM); % [MeV] % % - show Wmax -% ShowMe(Wmax_p,Ek_p,"W_{max} [MeV]","E_k [MeV]","W_{max} of PROTONs"); -% ShowMe(Wmax_C,Ek_C,"W_{max} [MeV]","E_k [MeV/u]","W_{max} of CARBON ions"); +% ShowMe(Wmax,myEk,"W_{max} [MeV]","E_k [MeV]",sprintf("W_{max} of %s",myPart)); % - density correction -densCorr_p_H2O=ComputeDensityCorrection(betagamma_p,densEff_x0_H2O,densEff_x1_H2O,densEff_a_H2O,densEff_m_H2O,densEff_C_H2O,densEff_d0_H2O); % [] -densCorr_C_H2O=ComputeDensityCorrection(betagamma_C,densEff_x0_H2O,densEff_x1_H2O,densEff_a_H2O,densEff_m_H2O,densEff_C_H2O,densEff_d0_H2O); % [] +densCorr_H2O=ComputeDensityCorrection(myBetaGamma,densEff_x0_H2O,densEff_x1_H2O,densEff_a_H2O,densEff_m_H2O,densEff_C_H2O,densEff_d0_H2O); % [] % - actual calculation -dEodx_p=ComputeBetheBloch(Zp,beta_p,betagamma_p,Wmax_p,ZoA_H2O,I_H2O,densCorr_p_H2O); % [MeV/g cm2] -dEodx_C=ComputeBetheBloch(ZC,beta_C,betagamma_C,Wmax_C,ZoA_H2O,I_H2O,densCorr_C_H2O); % [MeV/g cm2] +dEodx=ComputeBetheBloch(myZ,myBeta,myBetaGamma,Wmax,ZoA_H2O,I_H2O,densCorr_H2O); % [MeV/g cm2] % - show Bethe-Bloch -ShowMe(dEodx_p,betagamma_p," [MeV/g cm^2]","\beta\gamma []","Mean stopping power of PROTONs in WATER"); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); -ShowMe(dEodx_C,betagamma_C," [MeV/g cm^2]","\beta\gamma []","Mean stopping power of CARBON ions in WATER"); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); - -% - for selected energy: -% . relativistic quantities -[betaSel_p,gammaSel_p,betagammaSel_p]=ComputeRelativisticQuantities(EkSel_p,mp); % [], [], [] -[betaSel_C,gammaSel_C,betagammaSel_C]=ComputeRelativisticQuantities(EkSel_C*AC,mC); % [], [], [] -% . Wmax -WmaxSel_p=ComputeWmax(betagammaSel_p,gammaSel_p,mp); % [MeV] -WmaxSel_C=ComputeWmax(betagammaSel_C,gammaSel_C,mC); % [MeV] -% . density correction -densCorrSel_p_H2O=ComputeDensityCorrection(betagammaSel_p,densEff_x0_H2O,densEff_x1_H2O,densEff_a_H2O,densEff_m_H2O,densEff_C_H2O,densEff_d0_H2O); % [] -densCorrSel_C_H2O=ComputeDensityCorrection(betagammaSel_C,densEff_x0_H2O,densEff_x1_H2O,densEff_a_H2O,densEff_m_H2O,densEff_C_H2O,densEff_d0_H2O); % [] -% . actual calculation -dEodxSel_p=ComputeBetheBloch(Zp,betaSel_p,betagammaSel_p,WmaxSel_p,ZoA_H2O,I_H2O,densCorrSel_p_H2O); % [MeV/g cm2] -dEodxSel_C=ComputeBetheBloch(ZC,betaSel_C,betagammaSel_C,WmaxSel_C,ZoA_H2O,I_H2O,densCorrSel_C_H2O); % [MeV/g cm2] +ShowMe(dEodx,myBetaGamma," [MeV/g cm^2]","\beta\gamma []",sprintf("Mean stopping power of %s in WATER",myPart)); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); +ShowMe(dEodx,Ek," [MeV/g cm^2]",sprintf("E_k [%s]",unitEk),sprintf("Mean stopping power of %s in WATER",myPart)); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); %% compute range (based on Bethe-Bloch) -range_p=ComputeRange(Ek_p,dEodx_p*rho_H2O)*10; % [mm] -range_C=ComputeRange(Ek_C*AC,dEodx_C*rho_H2O)*10; % [mm] +range=ComputeRange(myEk,dEodx*rho_H2O)*10; % [mm] % - show range -ShowMe(range_p,Ek_p,"R [mm]","E_k [MeV]","Range of PROTONs in WATER"); -ShowMe(range_C,Ek_C,"R [mm]","E_k [MeV/u]","Range of CARBON ions in WATER"); -% - for selected energy: -rangeSel_p=interp1(betagamma_p,range_p,betagammaSel_p); % [mm] -rangeSel_C=interp1(betagamma_C,range_C,betagammaSel_C); % [mm] +ShowMe(range,Ek,"R [mm]",sprintf("E_k [%s]",unitEk),sprintf("Range of %s in WATER",myPart)); %% compute delta_p/p for n-mm of water-equivalent material (based on Bethe-Bloch) if (length(mmEquiv)==1) - % vs beam energy for a specific thickness % - actual calculation - dPoP_p=interp1(range_p,betagamma_p,range_p-mmEquiv)./betagamma_p-1; - dPoP_C=interp1(range_C,betagamma_C,range_C-mmEquiv)./betagamma_C-1; - % - show dPoP - ShowMe(dPoP_p,Ek_p,"\delta []","E_k [MeV]",sprintf("PROTON momentum variation after traversing %g mm of water equivalent",mmEquiv)); - ShowMe(dPoP_C,Ek_C,"\delta []","E_k [MeV/u]",sprintf("CARBON ion momentum variation after traversing %g mm of water equivalent",mmEquiv)); + dPoP=interp1(range,myBetaGamma,range-mmEquiv)./myBetaGamma-1; + % - show dPoP vs beam energy for a specific thickness + ShowMe(dPoP,Ek,"\delta []",sprintf("E_k [%s]",unitEk),sprintf("%s momentum variation after traversing %g mm of water equivalent",myPart,mmEquiv)); else - % vs thickness for a specific beam energy % - actual calculation - dPoP_p=interp1(range_p,betagamma_p,rangeSel_p-mmEquiv)./betagammaSel_p-1; - dPoP_C=interp1(range_C,betagamma_C,rangeSel_C-mmEquiv)./betagammaSel_C-1; - % - show dPoP - ShowMe(dPoP_p,mmEquiv,"\delta []","z_{H_2O} [mm]",sprintf("Momentum variation vs water equivalent thickness for %g MeV PROTONs",EkSel_p)); - ShowMe(dPoP_C,mmEquiv,"\delta []","z_{H_2O} [mm]",sprintf("Momentum variation vs water equivalent thickness for %g MeV/u CARBON IONs",EkSel_C)); + dPoP=interp1(range,myBetaGamma,range(end)-mmEquiv)./myBetaGamma(end)-1; + % - show dPoP vs thickness for a specific beam energy + ShowMe(dPoP,mmEquiv,"\delta []","z_{H_2O} [mm]",sprintf("Momentum variation vs water equivalent thickness for %g %s %s",Ek(end),unitEk,myPart)); end %% Landau-Vavilov for n-mm of water-equivalent material if (length(mmEquiv)==1) - % vs beam energy for a specific thickness % - actual calculation - mpEnLoss_p=ComputeLandauVavilov(Zp,beta_p,betagamma_p,mmEquiv/10*rho_H2O,ZoA_H2O,I_H2O,densCorr_p_H2O); % [MeV] - mpEnLoss_C=ComputeLandauVavilov(ZC,beta_C,betagamma_C,mmEquiv/10*rho_H2O,ZoA_H2O,I_H2O,densCorr_C_H2O); % [MeV] - % - show Landau-Vavilov - ShowMe(mpEnLoss_p/(mmEquiv/10*rho_H2O),betagamma_p,"\DeltaE [MeV/g cm^2]","\beta\gamma []",sprintf("Most probable energy loss of PROTONs after traversing %g mm of water equivalent",mmEquiv)); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); - ShowMe(mpEnLoss_C/(mmEquiv/10*rho_H2O),betagamma_C,"\DeltaE [MeV/g cm^2]","\beta\gamma []",sprintf("Most probable energy loss of CARBON IONs after traversing %g mm of water equivalent",mmEquiv)); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); + mpEnLoss=ComputeLandauVavilov(myZ,myBeta,myBetaGamma,mmEquiv/10*rho_H2O,ZoA_H2O,I_H2O,densCorr_H2O); % [MeV] + % - show Landau-Vavilov vs beam energy for a specific thickness + ShowMe(mpEnLoss/(mmEquiv/10*rho_H2O),myBetaGamma,"\DeltaE [MeV/g cm^2]","\beta\gamma []",sprintf("Most probable energy loss of %s after traversing %g mm of water equivalent",myPart,mmEquiv)); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); else - % vs thickness for a specific beam energy % - actual calculation - mpEnLoss_p=ComputeLandauVavilov(Zp,betaSel_p,betagammaSel_p,mmEquiv/10*rho_H2O,ZoA_H2O,I_H2O,densCorrSel_p_H2O); % [MeV] - mpEnLoss_C=ComputeLandauVavilov(ZC,betaSel_C,betagammaSel_C,mmEquiv/10*rho_H2O,ZoA_H2O,I_H2O,densCorrSel_C_H2O); % [MeV] - % - show Landau-Vavilov - ShowMe(mpEnLoss_p,mmEquiv,"\DeltaE [MeV]","z_{H_2O} [mm]",sprintf("Most probable energy loss vs water equivalent thickness for %g MeV PROTONs",EkSel_p)); set(gca, 'YScale', 'log'); % set(gca, 'XScale', 'log'); - ShowMe(mpEnLoss_C,mmEquiv,"\DeltaE [MeV]","z_{H_2O} [mm]",sprintf("Most probable energy loss vs water equivalent thickness for %g MeV/u CARBON IONs",EkSel_C)); set(gca, 'YScale', 'log'); % set(gca, 'XScale', 'log'); + mpEnLoss=ComputeLandauVavilov(myZ,myBeta(end),myBetaGamma(end),mmEquiv/10*rho_H2O,ZoA_H2O,I_H2O,densCorr_H2O(end)); % [MeV] + % - show Landau-Vavilov vs thickness for a specific beam energy + ShowMe(mpEnLoss,mmEquiv,"\DeltaE [MeV]","z_{H_2O} [mm]",sprintf("Most probable energy loss vs water equivalent thickness for %g %s %s IONs",Ek(end),unitEk,myPart)); set(gca, 'YScale', 'log'); % set(gca, 'XScale', 'log'); end %% compare Bethe-Bloch and Landau-Vavilov if (length(mmEquiv)==1) - ShowMyContent=NaN(2,size(dEodx_p,2)); ShowMyContent(1,:)=dEodx_p; ShowMyContent(2,:)=mpEnLoss_p/(mmEquiv/10*rho_H2O); myLegend=[" (Bethe-Bloch)" "most probable \DeltaE/x (Landau-Vavilov)"]; -% ShowMe(ShowMyContent,betagamma_p,"\DeltaE [MeV/g cm^2]","\beta\gamma []",sprintf("Bethe-Bloch vs Landau-Vavilov for PROTONs in %g mm of WATER",mmEquiv),myLegend); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); - ShowMe(ShowMyContent,Ek_p,"\DeltaE [MeV/g cm^2]","E_k [MeV]",sprintf("Bethe-Bloch vs Landau-Vavilov for PROTONs in %g mm of WATER",mmEquiv),myLegend); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); - ShowMyContent=NaN(2,size(dEodx_C,2)); ShowMyContent(1,:)=dEodx_C; ShowMyContent(2,:)=mpEnLoss_C/(mmEquiv/10*rho_H2O); myLegend=[" (Bethe-Bloch)" "most probable \DeltaE/x (Landau-Vavilov)"]; -% ShowMe(ShowMyContent,betagamma_C,"\DeltaE [MeV/g cm^2]","\beta\gamma []",sprintf("Bethe-Bloch vs Landau-Vavilov for CARBON ions in %g mm of WATER",mmEquiv),myLegend); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); - ShowMe(ShowMyContent,Ek_C,"\DeltaE [MeV/g cm^2]","E_k [MeV/u]",sprintf("Bethe-Bloch vs Landau-Vavilov for CARBON ions in %g mm of WATER",mmEquiv),myLegend); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); + ShowMyContent=NaN(2,size(dEodx,2)); ShowMyContent(1,:)=dEodx; ShowMyContent(2,:)=mpEnLoss/(mmEquiv/10*rho_H2O); myLegend=[" (Bethe-Bloch)" "most probable \DeltaE/x (Landau-Vavilov)"]; + ShowMe(ShowMyContent,myBetaGamma,"\DeltaE [MeV/g cm^2]","\beta\gamma []",sprintf("Bethe-Bloch vs Landau-Vavilov for %s in %g mm of WATER",myPart,mmEquiv),myLegend); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); + ShowMe(ShowMyContent,Ek,"\DeltaE [MeV/g cm^2]",sprintf("E_k [%s]",unitEk),sprintf("Bethe-Bloch vs Landau-Vavilov for %s in %g mm of WATER",myPart,mmEquiv),myLegend); set(gca, 'YScale', 'log'); set(gca, 'XScale', 'log'); else - dE_p=EkSel_p-interp1(range_p,Ek_p,rangeSel_p-mmEquiv); - ShowMyContent=NaN(2,size(dE_p,2)); ShowMyContent(1,:)=dE_p; ShowMyContent(2,:)=mpEnLoss_p; myLegend=[" (Bethe-Bloch)" "most probable \DeltaE (Landau-Vavilov)"]; - ShowMe(ShowMyContent,mmEquiv,"\DeltaE [MeV]","z_{H_2O} [mm]",sprintf("Bethe-Bloch vs Landau-Vavilov for %g MeV PROTONs in WATER",EkSel_p),myLegend); set(gca, 'YScale', 'log'); % set(gca, 'XScale', 'log'); - dE_C=EkSel_C-interp1(range_C,Ek_C,rangeSel_C-mmEquiv); dE_C=dE_C*AC; - ShowMyContent=NaN(2,size(dE_C,2)); ShowMyContent(1,:)=dE_C; ShowMyContent(2,:)=mpEnLoss_C; myLegend=[" (Bethe-Bloch)" "most probable \DeltaE (Landau-Vavilov)"]; - ShowMe(ShowMyContent,mmEquiv,"\DeltaE [MeV]","z_{H_2O} [mm]",sprintf("Bethe-Bloch vs Landau-Vavilov for %g MeV/u CARBON ions in WATER",EkSel_C),myLegend); set(gca, 'YScale', 'log'); % set(gca, 'XScale', 'log'); + dE=Ek(end)-interp1(range,Ek,range(end)-mmEquiv); + ShowMyContent=NaN(2,size(dE,2)); ShowMyContent(1,:)=dE; ShowMyContent(2,:)=mpEnLoss; myLegend=[" (Bethe-Bloch)" "most probable \DeltaE (Landau-Vavilov)"]; + ShowMe(ShowMyContent,mmEquiv,"\DeltaE [MeV]","z_{H_2O} [mm]",sprintf("Bethe-Bloch vs Landau-Vavilov for %g %s %s in WATER",Ek(end),unitEk,myPart),myLegend); set(gca, 'YScale', 'log'); % set(gca, 'XScale', 'log'); end %% local functions diff --git a/educational/ShowMCS.m b/educational/ShowMCS.m index eef4284..0f8814f 100644 --- a/educational/ShowMCS.m +++ b/educational/ShowMCS.m @@ -17,11 +17,12 @@ % - particle myPart="PROTON"; % available: "PROTON", "CARBON", "HELIUM" +% NOTA BENE: please choose either an array of values of Ek and a single value +% for the range traversed or the other way around; % - kinetic energies Ek=1:1:250; % [MeV] % proton energies % Ek=1:1:400; % [MeV/A] % carbon energies -% Ek_p=228.57; % single energy - +% Ek=228.57; % single energy % - range traversed mmEquiv=100.0; % [mm] % mmEquiv=0.1:0.1:10.1; % [mm] @@ -30,18 +31,8 @@ X0_H2O=36.08; % [cm] %% Load particle data - -run(".\particleData.m"); -switch upper(myPart) - case {"P","PROT","PROTON"} - myM=mp; myEk=Ek; myZ=Zp; unitEk="MeV"; - case {"HE","HELIUM","ALPHA","ALFA"} - myM=mHe; myEk=Ek*AHe; myZ=ZHe; unitEk="MeV/u"; - case {"C","CARB","CARBON"} - myM=mC; myEk=Ek*AC; myZ=ZC; unitEk="MeV/u"; - otherwise - error("Unknown particle %s!",myPart); -end +% returns: myM [MeV/c2], myEk [MeV], myZ [], unitEk ("MeV" for protons, "MeV/u" for others); +run(".\setParticle.m"); %% compute MCS scattering tables % - relativistic quantities @@ -58,13 +49,13 @@ end %% show effect of scattering -L=0.15; % distance travelled by scattered bea, [m] +L=0.15; % distance in vacuum travelled by scattered beam [m] if (length(mmEquiv)==1 && length(Ek)>1) % as a function of energy ShowMe(theta0*L*1000,Ek,"\theta\timesL [mm]",sprintf("E_k [%s]",unitEk),sprintf("MCS for %s after traversing %g mm of water equivalent",myPart,mmEquiv)); set(gca, 'YScale', 'log'); else % as a function of material thickness - ShowMe(theta0*1000,mmEquiv,"\theta\timesL [mm]","H_2O Range [mm]",sprintf("MCS for %s of %g %s",myPart,Ek_p,unitEk)); set(gca, 'YScale', 'log'); + ShowMe(theta0*1000,mmEquiv,"\theta\timesL [mm]","H_2O Range [mm]",sprintf("MCS for %s of %g %s",myPart,Ek,unitEk)); set(gca, 'YScale', 'log'); end diff --git a/educational/particleData.m b/educational/particleData.m index aabbd09..d3f3585 100644 --- a/educational/particleData.m +++ b/educational/particleData.m @@ -2,17 +2,24 @@ %% description % file with standard particle data +fprintf("reading standard particle data..."); + % proton data +fprintf("...PROTON data..."); Ap=1; Zp=1; mp=938.27208816; % [MeV/c2] % helium data +fprintf("...HELIUM_4^2+ data..."); AHe=4; ZHe=2; mHe=931.2386834*AHe; % [MeV/c2] % carbon data +fprintf("...CARBON_12^6+ data..."); AC=12; ZC=6; mC=931.2386834*AC; % [MeV/c2] + +fprintf("...done;"); diff --git a/educational/setParticle.m b/educational/setParticle.m new file mode 100644 index 0000000..c3184f9 --- /dev/null +++ b/educational/setParticle.m @@ -0,0 +1,17 @@ +%% description +% file for loading the correct particle data + +run(".\particleData.m"); +switch upper(myPart) + case {"P","PROT","PROTON"} + fprintf("...using PROTON data..."); + myM=mp; myEk=Ek; myZ=Zp; unitEk="MeV"; + case {"HE","HELIUM","ALPHA","ALFA"} + fprintf("...using HELIUM_4^2+ data..."); + myM=mHe; myEk=Ek*AHe; myZ=ZHe; unitEk="MeV/u"; + case {"C","CARB","CARBON"} + fprintf("...using CARBON_12^6+ data..."); + myM=mC; myEk=Ek*AC; myZ=ZC; unitEk="MeV/u"; + otherwise + error("Unknown particle %s!",myPart); +end From b622dfe7644418ee81bcad125f07e3aff5f25966 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Mon, 6 Nov 2023 14:41:35 +0100 Subject: [PATCH 15/23] creating a displayLib folder with functions used only by Display* scritps --- DisplayLGENs.m | 27 ++++++++------------------- DisplayRefValues.m | 14 +++++++------- README.md | 1 + displayLib/ConfigCheck.m | 30 ++++++++++++++++++++++++++++++ 4 files changed, 46 insertions(+), 26 deletions(-) create mode 100644 displayLib/ConfigCheck.m diff --git a/DisplayLGENs.m b/DisplayLGENs.m index 5a064b3..afd82cc 100644 --- a/DisplayLGENs.m +++ b/DisplayLGENs.m @@ -14,21 +14,22 @@ % ------------------------------------------------------------------------- % USER's input data -machine=["ISO1" "ISO2" "ISO3" "ISO4" ]; % "ISO1" "ISO2" "ISO3" "ISO4"]; +machine=["Sala2H" "Sala2H"];%["ISO1" "ISO2" "ISO3" "ISO4" ]; % "ISO1" "ISO2" "ISO3" "ISO4"]; +machine="ISO3";%["ISO1" "ISO2" "ISO3" "ISO4" ]; % "ISO1" "ISO2" "ISO3" "ISO4"]; % beamPart=["PROTON" "PROTON" "PROTON" "PROTON" ]; %"CARBON" "CARBON" "CARBON" "CARBON" ]; -beamPart=["CARBON" "CARBON" "CARBON" "CARBON" ]; +beamPart=["PROTON" "CARBON" ]; % machine="ISO2"; % beamPart=["PROTON" "CARBON" ]; config="TM"; % select configuration: TM, RFKO -myTitle="LGEN values as downloaded on 2023-01-30"; +myTitle="LGEN values"; % ------------------------------------------------------------------------- % ------------------------------------------------------------------------- % check of user input data nSets=max([length(beamPart) length(machine) length(config)]); -beamPart=MyCheck(beamPart,nSets,"beamPart"); -machine=MyCheck(machine,nSets,"machine"); -config=MyCheck(config,nSets,"config"); +beamPart=ConfigCheck(beamPart,nSets,"beamPart"); +machine=ConfigCheck(machine,nSets,"machine"); +config=ConfigCheck(config,nSets,"config"); %% clear variables [cyCodes,ranges,Eks,Brhos,currents,fields,kicks,psNames,FileNameCurrents,magNames]=... @@ -72,21 +73,9 @@ %% visual checks LGENvisualCheck(psNames,Eks ,"Ek [MeV/u]",currents,"I [A]",magNames,myTitle,myLeg); LGENvisualCheck(psNames,Eks ,"Ek [MeV/u]",normCurrents,"I/B\rho [A/Tm]",magNames,myTitle,myLeg); -LGENvisualCheck(psNames,ranges,"range [mm]",currents,"I [A]",magNames,myTitle,myLeg); +LGENvisualCheck(psNames,ranges,"range [mm]",currents,"I [A]",magNames,myTitle,myLeg,"QUE"); LGENvisualCheck(psNames,ranges,"range [mm]",normCurrents,"I/B\rho [A/Tm]",magNames,myTitle,myLeg); LGENvisualCheck(psNames,Brhos ,"B\rho [Tm]",currents,"I [A]",magNames,myTitle,myLeg); LGENvisualCheck(psNames,Brhos ,"B\rho [Tm]",normCurrents,"I/B\rho [A/Tm]",magNames,myTitle,myLeg); %% local functions -function OutVar=MyCheck(InVar,nSets,myName) - OutVar=InVar; - nIn=length(OutVar); - if (nIn1 and length(InVar)==1, +% the function replicates the single +% value of InVar nSets times; +% +% input +% - InVar (1D strings): a list of strings (eg beam particles, +% configurations, etc...); +% - nSets (scalar): expected length(InVar); +% - myName (string): a string used just for printout on terminal in case of +% error; +% +% output: +% - OutVar (1D strings): a list of strings such that length(OutVar)=nSets; +% + OutVar=InVar; + nIn=length(OutVar); + if (nIn Date: Mon, 6 Nov 2023 14:42:59 +0100 Subject: [PATCH 16/23] forgotten function in DisplayRefValues.m --- DisplayRefValues.m | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/DisplayRefValues.m b/DisplayRefValues.m index 1459465..4deddaa 100644 --- a/DisplayRefValues.m +++ b/DisplayRefValues.m @@ -125,16 +125,3 @@ function ShowStuff(xVals,yVals,showXlabel,showYlabel,myLeg,nCoupled,myTitle) end end - -function OutVar=ConfigCheck(InVar,nSets,myName) - OutVar=InVar; - nIn=length(OutVar); - if (nIn Date: Wed, 15 Nov 2023 16:44:33 +0100 Subject: [PATCH 17/23] bug fix concerning handling of time arrays when parsing beam profiles --- measurement_analysis/ParseBeamProfiles.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/measurement_analysis/ParseBeamProfiles.m b/measurement_analysis/ParseBeamProfiles.m index eddda8b..7e9da0d 100644 --- a/measurement_analysis/ParseBeamProfiles.m +++ b/measurement_analysis/ParseBeamProfiles.m @@ -300,7 +300,7 @@ end cyCodes=cyCodes(idx); cyProgs=string(cyProgs); - times=times(:,idx); + if (size(times,2)==nAcquired), times=times(:,idx); end else measData=missing; cyProgs=missing; From 30eb2e889f20f29db36eda3beb997db2118911b8 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 15 Nov 2023 16:47:30 +0100 Subject: [PATCH 18/23] Extracting some general-purpose functions from Display scripts --- DisplayRefValues.m | 68 +------------------------- general/ExtractFromTable.m | 21 ++++++++ general/ShowSeries.m | 99 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 122 insertions(+), 66 deletions(-) create mode 100644 general/ExtractFromTable.m create mode 100644 general/ShowSeries.m diff --git a/DisplayRefValues.m b/DisplayRefValues.m index 4deddaa..5f4fa99 100644 --- a/DisplayRefValues.m +++ b/DisplayRefValues.m @@ -50,7 +50,7 @@ FullFileName=ReturnDefFile("BRHO",myConfig); EnData=readtable(FullFileName); myLeg(iSet)=myConfig; % - store data - [tmpDataX,tmpDataY]=ExtractData(EnData,showX,showY); + [tmpDataX,tmpDataY]=ExtractFromTable(EnData,showX,showY); xVals=ExpandMat(xVals,tmpDataX); yVals=ExpandMat(yVals,tmpDataY); end @@ -58,70 +58,6 @@ yVals=permute(yVals,[1 3 2]); %% show stuff -ShowStuff(xVals,yVals,showXlabel,showYlabel,myLeg,nCoupled); +ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg); %% local functions -function ShowStuff(xVals,yVals,showXlabel,showYlabel,myLeg,nCoupled,myTitle) - if (~exist("nCoupled","var")), nCoupled=1; end - if (~exist("myTitle","var")), myTitle=missing(); end - nShows=size(xVals,3); - nSets=size(xVals,2); - - figure(); - cm=colormap(parula(nSets+1)); - % - set grid of plots - nPlots=nShows; - if (~ismissing(myLeg)), nPlots=nPlots+1; end - [nRows,nCols]=GetNrowsNcols(nPlots,nCoupled); - % - set markers - if ( nSets==2 ) - markers=[ "o" "*" ]; - else - markers=strings(nSets,1); - markers(:)="."; - end - % - actually plot - for iPlot=1:nShows - subplot(nRows,nCols,iPlot); - lFirst=true; - for iSet=1:nSets - if (lFirst), lFirst=false; else, hold on; end - plot(xVals(:,iSet,iPlot),yVals(:,iSet,iPlot),"-","Marker",markers(iSet),"Color",cm(iSet,:)); - end - xlabel(showXlabel(iPlot)); ylabel(showYlabel(iPlot)); grid on; - end - % - legend plot - if (~ismissing(myLeg)) - subplot(nRows,nCols,nPlots); - lFirst=true; - for iSet=1:nSets - if (lFirst), lFirst=false; else, hold on; end - plot(NaN(),NaN(),"-","Marker",markers(iSet),"Color",cm(iSet,:)); - end - legend(myLeg); - end - % - general stuff - if (~ismissing(myTitle)), sgtitle(myTitle); end -end - -function [DataX,DataY]=ExtractData(EnData,showX,showY) - names=fieldnames(EnData); - DataX=missing(); DataY=missing(); - for ii=1:length(showX) - isField=strcmpi(showX(ii),names); - if any(isField) - DataX=ExpandMat(DataX,EnData.(names{isField})); - else - error("no field %s in table!",showX(ii)); - end - end - for ii=1:length(showY) - isField=strcmpi(showY(ii),names); - if any(isField) - DataY=ExpandMat(DataY,EnData.(names{isField})); - else - error("no field %s in table!",showY(ii)); - end - end - -end diff --git a/general/ExtractFromTable.m b/general/ExtractFromTable.m new file mode 100644 index 0000000..8251927 --- /dev/null +++ b/general/ExtractFromTable.m @@ -0,0 +1,21 @@ +function [DataX,DataY]=ExtractFromTable(myTable,myXs,myYs) +% ExtractFromTable to get specific columns as X and Y values + fNames=fieldnames(myTable); + DataX=missing(); DataY=missing(); + for ii=1:length(myXs) + iField=strcmpi(myXs(ii),fNames); + if any(iField) + DataX=ExpandMat(DataX,myTable.(fNames{iField})); + else + error("no field %s in table!",myXs(ii)); + end + end + for ii=1:length(myYs) + iField=strcmpi(myYs(ii),fNames); + if any(iField) + DataY=ExpandMat(DataY,myTable.(fNames{iField})); + else + error("no field %s in table!",myYs(ii)); + end + end +end diff --git a/general/ShowSeries.m b/general/ShowSeries.m new file mode 100644 index 0000000..2082953 --- /dev/null +++ b/general/ShowSeries.m @@ -0,0 +1,99 @@ +function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,refLeg,nCoupled) +% ShowStuff to display data series (xVals,yVals) as several plots in a single figure +% +% NOTA BENE: +% - size(xVals)=size(yVals) (NOT CHECKED YET!) +% - number of points in each plotted series: size(xVals,1)=size(yVals,1); +% - number of series per plot: size(xVals,2)=size(yVals,2); +% - number of plots per figure: size(xVals,3)=size(yVals,3); +% - ref data can be: +% . unique for all plots or a dedicated one for each plot -- size(refY,3); +% . for each plot, a unique series (black) or a few (the first one black, the others grey) -- size(refY,2); + + if (~exist("myTitles","var")), myTitles=missing(); end + if (~exist("refY","var")), refY=missing(); end + if (~exist("refX","var")), refX=missing(); end + if (~exist("refLeg","var")), refLeg=missing(); end + if (~exist("nCoupled","var")), nCoupled=1; end + nSets=size(yVals,2); + nShows=size(yVals,3); + nPlotsLegendThresh=10; % max number of subplots before having a dedicated one for the legend + nSeriesLegendThresh=10; % max number of series before having a dedicated one for the legend + + figure(); + % - set grid of plots + nPlots=nShows; + if (all(~ismissing(myLeg)) && ( nPlots>nPlotsLegendThresh || nSets>nSeriesLegendThresh ) ), nPlots=nPlots+1; end + [nRows,nCols]=GetNrowsNcols(nPlots,nCoupled); + % - set markers + if ( nSets==2 ) + markers=[ "o" "*" ]; + cm=[1 0 0 ; 0 0 1]; % red and blue + else + markers=strings(nSets,1); markers(:)="."; + cm=colormap(parula(nSets)); + end + if (~ismissing(refY)) + refMarkers=strings(size(refY,2),1); refMarkers(:)="-"; + refColors=NaN(size(refY,2),3); refColors(1,:)=[0 0 0]; + if (size(refY,2)>1), refColors(2:end,:)=repmat([0.7 0.7 0.7],[size(refColors,1)-1 1]); end + end + % - legend + tmpLeg=[]; + if (all(~ismissing(myLeg))) + tmpLeg=myLeg; + end + if (all(~ismissing(refLeg))) + tmpLeg=[tmpLeg refLeg]; + end + % - actually plot + for iPlot=1:nShows + axs(iPlot)=subplot(nRows,nCols,iPlot); + lFirst=true; + if (~ismissing(refX)) + for jRef=1:size(refY,2) + if (lFirst), lFirst=false; else, hold on; end + if (size(refX,2)==1), jX=1; else, jX=jRef; end + if (size(refX,3)==1), kX=1; else, kX=iPlot; end + if (size(refY,2)==1), jY=1; else, jY=jRef; end + if (size(refY,3)==1), kY=1; else, kY=iPlot; end + plot(refX(:,jX,kX),refY(:,jY,kY),refMarkers(jRef),"Color",refColors(jRef,:)); + end + lFirst=false; + end + for iSet=1:nSets + if (lFirst), lFirst=false; else, hold on; end + if (size(xVals,2)==1), jX=1; else, jX=iSet; end + if (size(xVals,3)==1), kX=1; else, kX=iPlot; end + plot(xVals(:,jX,kX),yVals(:,iSet,iPlot),"-","Marker",markers(iSet),"Color",cm(iSet,:)); + end + if (length(showXlabel)==nShows), iX=iPlot; else, iX=1; end; xlabel(showXlabel(iX)); + if (length(showYlabel)==nShows), iY=iPlot; else, iY=1; end; ylabel(showYlabel(iY)); + grid on; + if (all(~ismissing(myTitles))), title(myTitles(iPlot)); end + if (~isempty(tmpLeg) && ( nPlots<=nPlotsLegendThresh && nSets<=nSeriesLegendThresh ) ), legend(tmpLeg); end + end + % - legend plot + if (~isempty(tmpLeg) && ( nPlots>nPlotsLegendThresh || nSets>nSeriesLegendThresh ) ) + subplot(nRows,nCols,nPlots); + lFirst=true; + if (all(~ismissing(refLeg))) + plot(NaN(),NaN(),"k-"); + if (size(refY,2)>1) + hold on; plot(NaN(),NaN(),"-","Color",[.7 .7 .7]); + hold on; plot(NaN(),NaN(),"-","Color",[.7 .7 .7]); + end + lFirst=false; + end + for iSet=1:nSets + if (lFirst), lFirst=false; else, hold on; end + plot(NaN(),NaN(),"-","Marker",markers(iSet),"Color",cm(iSet,:)); + end + legend(tmpLeg,"Location","best"); title("legend"); + end + % - general stuff + if (all(~ismissing(myTitles))) + if (length(myTitles)>nShows), sgtitle(myTitles(end)); end + end + % linkaxes(axs,"xy"); +end From 52ff55d716c5b1675d2b1e3f2b66b411296b32e8 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 15 Nov 2023 16:49:23 +0100 Subject: [PATCH 19/23] implementing handling of reference spot sizes --- DisplaySpots.m | 100 ++++++++++++++++++ measurement_analysis/ParseMPrefFWHMfile.m | 31 ++++++ measurement_analysis/Spots_FWHMRefSeries.m | 30 ++++++ .../Spots_LoadEffectiveSpecs.m | 31 ++++++ measurement_analysis/Spots_LoadRef.m | 21 ++++ measurement_analysis/Spots_MeritFWHM.m | 17 +++ operations/cyCodes/MapCyCodes.m | 26 +++++ 7 files changed, 256 insertions(+) create mode 100644 DisplaySpots.m create mode 100644 measurement_analysis/ParseMPrefFWHMfile.m create mode 100644 measurement_analysis/Spots_FWHMRefSeries.m create mode 100644 measurement_analysis/Spots_LoadEffectiveSpecs.m create mode 100644 measurement_analysis/Spots_LoadRef.m create mode 100644 measurement_analysis/Spots_MeritFWHM.m create mode 100644 operations/cyCodes/MapCyCodes.m diff --git a/DisplaySpots.m b/DisplaySpots.m new file mode 100644 index 0000000..59cbf49 --- /dev/null +++ b/DisplaySpots.m @@ -0,0 +1,100 @@ +% {}~ + +%% description +% this is a script which displays beam spots; + +%% include libraries +% - include Matlab libraries +pathToLibrary="./"; +addpath(genpath(pathToLibrary)); +pathToLibrary="../MachineRefs"; +addpath(genpath(pathToLibrary)); + +%% clear +clear all; +close all; + +%% settings + +% ------------------------------------------------------------------------- +% USER's input data +machine=["Sala1" "Sala2H" "Sala2V" "Sala3"]; +beamPart="CARBON"; +config="TM"; % select configuration: TM, RFKO +myXwhat="Range"; % show stuff as a function of +showXlabel="R_{H_2O} [mm]"; +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% check of user input data +nSets=max([length(beamPart) length(machine) length(config)]); +beamPart=ConfigCheck(beamPart,nSets,"beamPart"); +machine=ConfigCheck(machine,nSets,"machine"); +config=ConfigCheck(config,nSets,"config"); + +%% clear variables +[SP_Eks SP_Mms SP_FWHMs]=... + deal(missing(),missing(),missing()); +[MP_Eks MP_Mms MP_FWHMs MP_myPos]=... + deal(missing(),missing(),missing(),missing()); + +%% get values + +% - acquire Steering Pazienti stuff +for iSet=1:nSets + % - parse DB file + [tmpEks,tmpMms,tmpFWHMs]=Spots_LoadRef(machine(iSet),beamPart(iSet),config(iSet),"SteerPaz"); + % - store data + SP_Eks=ExpandMat(SP_Eks,repmat(tmpEks,[1 size(tmpFWHMs,2)])); + SP_Mms=ExpandMat(SP_Mms,repmat(tmpMms,[1 size(tmpFWHMs,2)])); + SP_FWHMs=ExpandMat(SP_FWHMs,tmpFWHMs); +end + +% - acquire Medical Physics stuff +for iSet=1:nSets + % - parse DB file + [tmpEks,tmpMms,tmpFWHMs,tmpMyPos]=Spots_LoadRef(machine(iSet),beamPart(iSet),config(iSet),"MedPhys","ALL"); + % - store data + MP_Eks=ExpandMat(MP_Eks,repmat(tmpEks,[1 size(tmpFWHMs,2)])); + MP_Mms=ExpandMat(MP_Mms,repmat(tmpMms,[1 size(tmpFWHMs,2)])); + MP_FWHMs=ExpandMat(MP_FWHMs,tmpFWHMs); + MP_myPos=ExpandMat(MP_myPos,tmpMyPos'); +end + +if (size(MP_FWHMs,3)~=size(SP_FWHMs,3)) + error("something wrong at parsing: %d MedPhys curves, and %d SteerPaz curves",size(MP_FWHMs,2),size(SP_FWHMs,3)); +end + +%% show data +nDataSets=length(machine); +% . labels for longitudinal position +PosLabs=strings(size(MP_myPos,1),1); +iIndices=(MP_myPos(:,1)-0.713<0); PosLabs(iIndices)=MP_myPos(iIndices,1)-0.713; PosLabs(iIndices)="ISO " +PosLabs(iIndices)+"m"; +iIndices=(MP_myPos(:,1)-0.713>0); PosLabs(iIndices)=MP_myPos(iIndices,1)-0.713; PosLabs(iIndices)="ISO +"+PosLabs(iIndices)+"m"; +iIndices=(MP_myPos(:,1)-0.713==0); PosLabs(iIndices)="ISO"; +% . labels for lines +LineLabs=strings(nDataSets,1); +LineLabs=machine(:)+" - "+beamPart(:); + +% - Steering Pazienti vs Medical Physics +myTitles=strings(nDataSets,1); myTitles(:)=LineLabs; myTitles(end+1)="SteerPaz vs MedPhys"; +iISO=(MP_myPos(:,1)==713E-3); +xRef=MP_Mms(:,iISO,:); yRef=MP_FWHMs(:,iISO,:); yRef=Spots_FWHMRefSeries(yRef); +ShowSeries(SP_Mms,SP_FWHMs,"R_{H_2O} [mm]","FWHM [mm]",["FWHM_x" "FWHM_y"],myTitles,yRef,xRef,["MP" "MP+" "MP-"]); + +% - Medical Physics: various positions (grouped by line) +myTitles=strings(nDataSets,1); myTitles(:)=LineLabs; myTitles(end)="MedPhys - various longitudinal positions"; +myLeg=PosLabs; +ShowSeries(MP_Mms,MP_FWHMs,"R_{H_2O} [mm]","FWHM [mm]",myLeg,myTitles); + +% - Medical Physics: various positions (grouped by position) +myTitles=strings(size(MP_myPos,1),1); myTitles(:)=PosLabs; myTitles(end+1)="MedPhys - various longitudinal positions"; +myLeg=LineLabs; +ShowSeries(permute(MP_Mms,[1 3 2]),permute(MP_FWHMs,[1 3 2]),"R_{H_2O} [mm]","FWHM [mm]",myLeg,myTitles); + +% - Medical Physics: various positions (grouped by position) +myTitles=strings(nDataSets,1); myTitles(:)=LineLabs; myTitles(end+1)="MedPhys - various longitudinal positions"; +myLeg=missing(); +ShowSeries(MP_myPos(:,1),permute(MP_FWHMs,[2 1 3]),"S [m]","FWHM [mm]",myLeg,myTitles); + +%% local functions diff --git a/measurement_analysis/ParseMPrefFWHMfile.m b/measurement_analysis/ParseMPrefFWHMfile.m new file mode 100644 index 0000000..1e15e19 --- /dev/null +++ b/measurement_analysis/ParseMPrefFWHMfile.m @@ -0,0 +1,31 @@ +function [Eks,Mms,FWHMs,myPos]=ParseMPrefFWHMfile(fileName,beamPart,loc) +% output vars: +% - size(FWHMs,1)=length(Ek)=lenght(Mms); +% - size(FWHMs,2)=length(myPos); + if (~exist("beamPart","var")), beamPart=missing(); end + if (~exist("loc","var")), loc=missing(); end + if (ismissing(beamPart)), beamPart="Prot"; end % default + if (ismissing(loc)), loc="ISO"; end % default + fprintf("parsing MedPhys ref file for FWHM values %s ...\n",fileName); + myData=readmatrix(fileName,"NumHeaderLines",1); + switch upper(beamPart) + case {"P","PROT","PROTON"} + Eks=myData(:,3); + case {"C","CARB","CARBON"} + Eks=myData(:,2); + otherwise + error("Unknown particle %s!",beamPart); + end + Mms=myData(:,4); + switch upper(loc) + case "ISO" + FWHMs=myData(:,8); + myPos=713E-3; % [m] + case "ALL" + FWHMs=myData(:,6:9); + myPos=[40 363 713 1213]*1E-3; % [mm] + otherwise + error("Unknown location to filter: %s",loc); + end + fprintf("...load data at %d locations;\n",length(myPos)); +end diff --git a/measurement_analysis/Spots_FWHMRefSeries.m b/measurement_analysis/Spots_FWHMRefSeries.m new file mode 100644 index 0000000..16519a2 --- /dev/null +++ b/measurement_analysis/Spots_FWHMRefSeries.m @@ -0,0 +1,30 @@ +function FWHMout=Spots_FWHMRefSeries(FWHMin,relLev,absLev) +% Spots_FWHMRefSeries to build numerical series representing +% reference values of FWHM, including +% tolerances (most stringent among relative +% and absolute ones); +% +% input vars: +% - FWHMin (float([nPoints,nSeries,nCases])): FWHM values [mm] (geo mean); +% - relLev (float): relative tolerance on reference values [0:1]; +% - absLev (float): absolute tolerance on reference values [mm; +% +% output vars: +% - FWHMout (float([nPoints,3,nCases])): ref FWHM values with tolerances [mm]: +% . FWHMout(:,1,:): actual reference values; +% . FWHMout(:,2,:): reference values + tolerance; +% . FWHMout(:,3,:): reference values - tolerance; +% +% + if (~exist("relLev","var")), relLev=0.1; end % default: 10% + if (~exist("absLev","var")), absLev=1.0; end % default: 1 mm + % - reference FWHM + if (size(FWHMin,2)==1) + FWHMout=FWHMin; + else + FWHMout=mean(FWHMin,2); + end + % - tolerance around reference + FWHMout(:,2,:)=min(FWHMout(:,1,:)*(1+relLev),FWHMout(:,1,:)+absLev); + FWHMout(:,3,:)=max(FWHMout(:,1,:)*(1-relLev),FWHMout(:,1,:)-absLev); +end diff --git a/measurement_analysis/Spots_LoadEffectiveSpecs.m b/measurement_analysis/Spots_LoadEffectiveSpecs.m new file mode 100644 index 0000000..2b97a00 --- /dev/null +++ b/measurement_analysis/Spots_LoadEffectiveSpecs.m @@ -0,0 +1,31 @@ +function [refFWHM,xVals]=Spots_LoadEffectiveSpecs(beamPart,myXwhat,config,spotType,relLev,absLev) + if (~exist("myXwhat","var")), myXwhat="Range"; end + if (~exist("config","var")), config="TM"; end + if (~exist("spotType","var")), spotType="STEERPAZ"; end + if (~exist("relLev","var")), relLev=0.1; end % default: 10% + if (~exist("absLev","var")), absLev=1.0; end % default: 1 mm + machine=["Sala1" "Sala2H" "Sala2V" "Sala3"]; + nSets=length(machine); + FWHMs=missing(); + for iSet=1:nSets + % - parse DB file + clear tmpCyProgs tmpCyCodes tmpBARs tmpFWHMs tmpASYMs tmpINTs tmpXvals; + myConfig=sprintf("%s,%s,%s",machine(iSet),beamPart,config); + FullFileName=ReturnDefFile(spotType,myConfig); + switch upper(spotType) + case "STEERPAZ" + [tmpCyProgs,tmpCyCodes,tmpBARs,tmpFWHMs,tmpASYMs,tmpINTs]=ParseBeamProfileSummaryFiles(FullFileName,"CAM"); + if (length(tmpCyProgs)<=1), error("...no ary data aquired!"); end + otherwise + error("Unknown spot type %s!",spotType); + end + % - convert cyCodes + [tmpXvals,~]=MapCyCodes(tmpCyCodes,myXwhat); + % - store data + FWHMs=ExpandMat(FWHMs,tmpFWHMs); + if (iSet==1), xVals=tmpXvals; end + end + % - compute actual reference + [FWHMgeoMean,refASYM]=Spots_MeritFWHM(FWHMs); + refFWHM=Spots_FWHMRefSeries(FWHMgeoMean,relLev,absLev); +end diff --git a/measurement_analysis/Spots_LoadRef.m b/measurement_analysis/Spots_LoadRef.m new file mode 100644 index 0000000..34c013c --- /dev/null +++ b/measurement_analysis/Spots_LoadRef.m @@ -0,0 +1,21 @@ +function [Eks,Mms,FWHMs,pos]=Spots_LoadRef(machine,beamPart,config,spotType,loc) + if (~exist("config","var")), config="TM"; end + if (~exist("spotType","var")), spotType="STEERPAZ"; end + if (~exist("loc","var")), loc="ISO"; end + myConfig=sprintf("%s,%s,%s",machine,beamPart,config); + FullFileName=ReturnDefFile(spotType,myConfig); + switch upper(spotType) + case {"STEERINGPAZIENTI","SP","STEERPAZ"} + [cyProgs,cyCodes,BARs,FWHMs,ASYMs,INTs]=ParseBeamProfileSummaryFiles(FullFileName,"CAM"); + if (length(cyProgs)<=1), error("...no data aquired!"); end + % - convert cyCodes + [Eks,~]=MapCyCodes(cyCodes,"Ek"); + [Mms,~]=MapCyCodes(cyCodes,"Range"); + pos=NaN(); + case {"MEDICALPHYSICS","MP","MEDPHYS"} + [Eks,Mms,FWHMs,pos]=ParseMPrefFWHMfile(FullFileName,beamPart,loc); + if (length(Eks)<=1), error("...no data aquired!"); end + otherwise + error("Unknown spot type %s!",spotType); + end +end diff --git a/measurement_analysis/Spots_MeritFWHM.m b/measurement_analysis/Spots_MeritFWHM.m new file mode 100644 index 0000000..a0dee70 --- /dev/null +++ b/measurement_analysis/Spots_MeritFWHM.m @@ -0,0 +1,17 @@ +function [FWHMgeoMean,ASYMs]=Spots_MeritFWHM(FWHMin) +% Spots_MeritFWHM to calculate some figures of merit of spot sizes +% +% input vars: +% - FWHMin (float([nPoints,nPlanes,nSeries])): FWHM values [mm]; +% +% output vars: +% - FWHMout (float([nPoints,nSeries])): geo mean FWHM values over planes [mm]; +% - ASYMs (float([nPoints,nSeries])): FWHM_y-FWHM_x values [mm]; +% + % - geo mean FWHM + FWHMgeoMean=sqrt(FWHMin(:,1,:).*FWHMin(:,2,:)); + FWHMgeoMean=reshape(FWHMgeoMean,[size(FWHMgeoMean,1),size(FWHMgeoMean,3)]); + % - xy-asymmetry + ASYMs=FWHMin(:,2,:)-FWHMin(:,1,:); + ASYMs=reshape(ASYMs,[size(ASYMs,1),size(ASYMs,3)]); +end diff --git a/operations/cyCodes/MapCyCodes.m b/operations/cyCodes/MapCyCodes.m new file mode 100644 index 0000000..5efc56e --- /dev/null +++ b/operations/cyCodes/MapCyCodes.m @@ -0,0 +1,26 @@ +function [vOut,Refs]=MapCyCodes(cyCodesIN,what,machine) + if (~exist("what","var")), what=missing(); end + if (~exist("machine","var")), machine=missing(); end + + if (ismissing(what)), what="Ek"; end + if (ismissing(machine)), machine="SYNCHRO"; end + + % - values ref values + [PcyCodes,whatValsP]=ExtractFromTable(readtable(ReturnDefFile("BRHO",sprintf("%s,PROTON",machine))),"CyCode",what); + [CcyCodes,whatValsC]=ExtractFromTable(readtable(ReturnDefFile("BRHO",sprintf("%s,CARBON",machine))),"CyCode",what); + + % - get mapping + [rangeCodes,partCodes]=DecodeCyCodes(cyCodesIN); + [lCCP,iCCP]=MapCyCode(rangeCodes,upper(PcyCodes)); + [lCCC,iCCC]=MapCyCode(rangeCodes,upper(CcyCodes)); + + % - assign values + vOut=NaN(length(cyCodesIN),1); + indicesP=(lCCP & FlagPart(partCodes,"p") ); + indicesC=(lCCC & FlagPart(partCodes,"C") ); + Refs=NaN(max(length(PcyCodes),length(CcyCodes)),2); + vOut(indicesP)=whatValsP(iCCP(indicesP)); + vOut(indicesC)=whatValsC(iCCC(indicesC)); + Refs(1:length(PcyCodes),1)=whatValsP; + Refs(1:length(CcyCodes),2)=whatValsC; +end From 0fe76e3611cd7934b5df8606d6f8b9ba20990390 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 15 Nov 2023 16:58:50 +0100 Subject: [PATCH 20/23] forgotten commit --- .../Spots_LoadEffectiveSpecs.m | 35 +++++++++++-------- 1 file changed, 20 insertions(+), 15 deletions(-) diff --git a/measurement_analysis/Spots_LoadEffectiveSpecs.m b/measurement_analysis/Spots_LoadEffectiveSpecs.m index 2b97a00..f181a53 100644 --- a/measurement_analysis/Spots_LoadEffectiveSpecs.m +++ b/measurement_analysis/Spots_LoadEffectiveSpecs.m @@ -9,23 +9,28 @@ FWHMs=missing(); for iSet=1:nSets % - parse DB file - clear tmpCyProgs tmpCyCodes tmpBARs tmpFWHMs tmpASYMs tmpINTs tmpXvals; - myConfig=sprintf("%s,%s,%s",machine(iSet),beamPart,config); - FullFileName=ReturnDefFile(spotType,myConfig); - switch upper(spotType) - case "STEERPAZ" - [tmpCyProgs,tmpCyCodes,tmpBARs,tmpFWHMs,tmpASYMs,tmpINTs]=ParseBeamProfileSummaryFiles(FullFileName,"CAM"); - if (length(tmpCyProgs)<=1), error("...no ary data aquired!"); end - otherwise - error("Unknown spot type %s!",spotType); - end - % - convert cyCodes - [tmpXvals,~]=MapCyCodes(tmpCyCodes,myXwhat); + [tmpEks,tmpMms,tmpFWHMs]=Spots_LoadRef(machine(iSet),beamPart,config,spotType); % - store data FWHMs=ExpandMat(FWHMs,tmpFWHMs); - if (iSet==1), xVals=tmpXvals; end + if (iSet==1) + switch upper(myXwhat) + case {"RANGE","R","MM"} + xVals=tmpMms; + case {"ENERGY","EK","MEV"} + xVals=tmpEks; + otherwise + error("...unable to understand request %s!",myXwhat); + end + end end % - compute actual reference - [FWHMgeoMean,refASYM]=Spots_MeritFWHM(FWHMs); - refFWHM=Spots_FWHMRefSeries(FWHMgeoMean,relLev,absLev); + switch upper(spotType) + case {"STEERINGPAZIENTI","SP","STEERPAZ"} + [refFWHM,refASYM]=Spots_MeritFWHM(FWHMs); + case {"MEDICALPHYSICS","MP","MEDPHYS"} + refFWHM=FWHMs; + otherwise + error("Unknown spot type %s!",spotType); + end + refFWHM=Spots_FWHMRefSeries(refFWHM,relLev,absLev); end From a8fc73e9a17a67366a7194708c2e4dec1803097d Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Wed, 15 Nov 2023 17:18:08 +0100 Subject: [PATCH 21/23] showing reference spot sizes when crunching CAM data --- DisplayBeamProfiles.m | 79 ++++++++++++++++++++++++++++++------------- general/ShowSeries.m | 8 ++--- 2 files changed, 59 insertions(+), 28 deletions(-) diff --git a/DisplayBeamProfiles.m b/DisplayBeamProfiles.m index 0bbe2a2..fc49b23 100644 --- a/DisplayBeamProfiles.m +++ b/DisplayBeamProfiles.m @@ -16,6 +16,10 @@ % for CAMeretta/DDS/GIM only, the script also compares summary data % against statistics data computed on profiles; +%% clean +clear all; +close all; + %% manual run if (~exist("MonPaths","var")) % script manually run by user @@ -27,6 +31,8 @@ if (~exist("pathToLibrary","var")) pathToLibrary=".\"; addpath(genpath(pathToLibrary)); + pathToLibrary="../MachineRefs"; + addpath(genpath(pathToLibrary)); end % - clear settings clear kPath myTit monTypes MonPaths myLabels @@ -35,33 +41,25 @@ % USER's input data % ------------------------------------------------------------------------- kPath="P:\Accelerating-System\Accelerator-data"; - monTypes="CAMdumps"; % CAM/CAMdumps, DDS, GIM, QBM/QPP/PIB/PMM/SFH/SFM/SFP + monTypes="CAM"; % CAM/CAMdumps, DDS, GIM, QBM/QPP/PIB/PMM/SFH/SFM/SFP myLabels=[... - "test 1: H scan (1E6 per spot)" - "test 2: grid (1E6 per spot)" - "test 3: V scan (1E6 per spot)" - "test 4: V scan (1.2E6 per spot)" - "test 5: V scan (1.2E6 per spot)" - "test 6: H scan (1.2E6 per spot)" - "test 7: H scan (1.2E6 per spot)" - "test 8: grid (1.2E6 per spot)" + "Sala1 (16-10-2023)" + "Sala2H (16-10-2023)" + "Sala2V (16-10-2023)" + "Sala3 (16-10-2023)" ]; % myLabels=monTypes; lSkip=false; % DDS summary file: skip first 2 lines (in addition to header line) myFigPath="."; % part-dependent stuff % - protoni - myFigName="Tests with DDS"; - myTit="Tests with DDS"; + myFigName="Machie Photos"; + myTit="Machine Photos"; MonPaths=[... - "P:\Accelerating-System\Accelerator-data\scambio\Alessio\2023-08-21_testsOcchiConiglio\DumpProtSO1_LineT_Size10_22-08-2023_2213\" - "P:\Accelerating-System\Accelerator-data\scambio\Alessio\2023-08-21_testsOcchiConiglio\DumpProtSO1_LineT_Size10_22-08-2023_2219\" - "P:\Accelerating-System\Accelerator-data\scambio\Alessio\2023-08-21_testsOcchiConiglio\DumpProtSO1_LineT_Size10_22-08-2023_2222\" - "P:\Accelerating-System\Accelerator-data\scambio\Alessio\2023-08-21_testsOcchiConiglio\DumpProtSO1_LineT_Size10_22-08-2023_2225\" - "P:\Accelerating-System\Accelerator-data\scambio\Alessio\2023-08-21_testsOcchiConiglio\DumpProtSO1_LineT_Size10_22-08-2023_2227\" - "P:\Accelerating-System\Accelerator-data\scambio\Alessio\2023-08-21_testsOcchiConiglio\DumpProtSO1_LineT_Size10_22-08-2023_2229\" - "P:\Accelerating-System\Accelerator-data\scambio\Alessio\2023-08-21_testsOcchiConiglio\DumpProtSO1_LineT_Size10_22-08-2023_2231\" - "P:\Accelerating-System\Accelerator-data\scambio\Alessio\2023-08-21_testsOcchiConiglio\DumpProtSO1_LineT_Size10_22-08-2023_2233\" + "P:\Accelerating-System\Accelerator-data\Area dati MD\00Steering\SteeringPazienti\carbonio\sala1\2023\2023.10.16\CarbSO2_LineZ_Size6_16-10-2023_0800\" + "P:\Accelerating-System\Accelerator-data\Area dati MD\00Steering\SteeringPazienti\carbonio\sala2U\2023\2023.10.16\CarbSO2_LineU_Size6_16-10-2023_0219\" + "P:\Accelerating-System\Accelerator-data\Area dati MD\00Steering\SteeringPazienti\carbonio\sala2V\2023\2023.10.16\CarbSO2_LineV_Size6_16-10-2023_0908\" + "P:\Accelerating-System\Accelerator-data\Area dati MD\00Steering\SteeringPazienti\carbonio\sala3\2023\2023.10.15\CarbSO2_LineT_Size6_15-10-2023_2313\" ]; % % - carbonio % myFigName="summary_carbonio_GIM_2023-05-09.10"; @@ -69,7 +67,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="ID"; % ["Ek"/"En"/"Energy","mm"/"r"/"range","ID"/"IDs"] + vsX="Ek"; % ["Ek"/"En"/"Energy","mm"/"r"/"range","ID"/"IDs"] iNotShow=false(127,2); iNotShow(1:2,1)=true; % do not show left-most fibers on hor plane (broken) end @@ -137,8 +135,8 @@ [tmpBARsProf,tmpFWHMsProf,tmpINTsProf]=StatDistributionsBDProcedure(tmpProfiles); end % - Eks,mms - tmpEksProf=ConvertCyCodes(tmpCyCodesProf,"Ek","MeVvsCyCo_P.xlsx"); - tmpMmsProf=ConvertCyCodes(tmpCyCodesProf,"mm","MeVvsCyCo_P.xlsx"); + tmpEksProf=MapCyCodes(tmpCyCodesProf,"Ek","SYNCHRO"); + tmpMmsProf=MapCyCodes(tmpCyCodesProf,"Range","SYNCHRO"); % - store data cyProgsProf=ExpandMat(cyProgsProf,tmpCyProgsProf); cyCodesProf=ExpandMat(cyCodesProf,tmpCyCodesProf); @@ -157,8 +155,8 @@ % - quick check of consistency of parsed data if (length(tmpCyProgsSumm)~=length(tmpCyProgsProf)), error("...inconsistent data set between summary data and actual profiles"); end % - Eks,mms - tmpEksSumm=ConvertCyCodes(tmpCyCodesSumm,"Ek","MeVvsCyCo_P.xlsx"); - tmpMmsSumm=ConvertCyCodes(tmpCyCodesSumm,"mm","MeVvsCyCo_P.xlsx"); + tmpEksSumm=MapCyCodes(tmpCyCodesSumm,"Ek","SYNCHRO"); + tmpMmsSumm=MapCyCodes(tmpCyCodesSumm,"Range","SYNCHRO"); % - store data cyProgsSumm=ExpandMat(cyProgsSumm,tmpCyProgsSumm); cyCodesSumm=ExpandMat(cyCodesSumm,tmpCyCodesSumm); @@ -222,6 +220,39 @@ end end +%% show figures of merit +% - references +[refFWHM,refXVals]=Spots_LoadEffectiveSpecs("CARBON",vsX,"TM"); refLeg=["reference" "ref+" "ref-"]; + +% - FWHM +showYlabels=strings(nDataSets,1); showYlabels(:)="FWHM [mm]"; +showXlabels=strings(nDataSets,1); showXlabels(:)=addLabel; +myTitles=strings(nDataSets+1,1); myTitles(1:end-1)=myLabels; myTitles(end)="FWHM"; +xVals=NaN(size(addIndex,1),2,size(addIndex,2)); +xVals(:,1,:)=addIndex; xVals(:,2,:)=addIndex; +myLeg=["HOR" "VER"]; +ShowSeries(xVals,FWHMsProf,showXlabels,showYlabels,myLeg,myTitles,refFWHM,refXVals,refLeg); + +% - xy-asymmetry +[FWHMprofGeoMean,profASYM]=Spots_MeritFWHM(FWHMsProf); +showYlabels=["[mm]" "[%]"]; +showXlabels=strings(2,1); showXlabels(:)=addLabel; +myLeg=myLabels; myTitles=["FWHM_y-FWHM_x" "normalised: (FWHM_y-FWHM_x)/geoAve" "yx-asymmetry"]; +xVals=NaN(size(addIndex,1),size(addIndex,2),2); +xVals(:,:,1)=addIndex; xVals(:,:,2)=addIndex; +yVals=NaN(size(addIndex,1),size(addIndex,2),2); +yVals(:,:,1)=profASYM; yVals(:,:,2)=profASYM./FWHMprofGeoMean*100; +ShowSeries(xVals,yVals,showXlabels,showYlabels,myLeg,myTitles); + +% - BARicenters +showYlabels=strings(nDataSets,1); showYlabels(:)="BAR [mm]"; +showXlabels=strings(nDataSets,1); showXlabels(:)=addLabel; +myTitles=strings(nDataSets+1,1); myTitles(1:end-1)=myLabels; myTitles(end)="BAR"; +xVals=NaN(size(addIndex,1),2,size(addIndex,2)); +xVals(:,1,:)=addIndex; xVals(:,2,:)=addIndex; +myLeg=["HOR" "VER"]; +ShowSeries(xVals,BARsProf,showXlabels,showYlabels,myLeg,myTitles); + %% save summary data % oFileName=strcat(kPath,"\scambio\Alessio\Carbonio_preSteering_summary-from-profiles.csv"); % SaveBeamProfileSummaryFile(oFileName,tmpBARsProf,tmpFWHMsProf,tmpINTsProf,tmpCyCodesProf,tmpCyProgsProf,"DDS"); diff --git a/general/ShowSeries.m b/general/ShowSeries.m index 2082953..c3edf64 100644 --- a/general/ShowSeries.m +++ b/general/ShowSeries.m @@ -40,11 +40,11 @@ function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,r end % - legend tmpLeg=[]; - if (all(~ismissing(myLeg))) - tmpLeg=myLeg; - end if (all(~ismissing(refLeg))) - tmpLeg=[tmpLeg refLeg]; + tmpLeg=refLeg; + end + if (all(~ismissing(myLeg))) + tmpLeg=[tmpLeg myLeg]; end % - actually plot for iPlot=1:nShows From 0568a585e07a67336aa470429b3179a8bba58ac9 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Thu, 16 Nov 2023 12:51:41 +0100 Subject: [PATCH 22/23] path to libraries added only if variable is not defined --- DisplaySpots.m | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/DisplaySpots.m b/DisplaySpots.m index 59cbf49..2925e3d 100644 --- a/DisplaySpots.m +++ b/DisplaySpots.m @@ -5,10 +5,12 @@ %% include libraries % - include Matlab libraries -pathToLibrary="./"; -addpath(genpath(pathToLibrary)); -pathToLibrary="../MachineRefs"; -addpath(genpath(pathToLibrary)); +if (~exist("pathToLibrary","var")) + pathToLibrary="./"; + addpath(genpath(pathToLibrary)); + pathToLibrary="../MachineRefs"; + addpath(genpath(pathToLibrary)); +end %% clear clear all; @@ -19,7 +21,7 @@ % ------------------------------------------------------------------------- % USER's input data machine=["Sala1" "Sala2H" "Sala2V" "Sala3"]; -beamPart="CARBON"; +beamPart="PROTON"; config="TM"; % select configuration: TM, RFKO myXwhat="Range"; % show stuff as a function of showXlabel="R_{H_2O} [mm]"; From 61e9415caf359066ff76c8c93af33276fcbdc6c2 Mon Sep 17 00:00:00 2001 From: Mereghetti Alessio Date: Mon, 4 Dec 2023 11:56:50 +0100 Subject: [PATCH 23/23] adding possibility to load sphinx data --- DisplaySpots.m | 29 +++++++-- general/ShowSeries.m | 60 +++++++++++++++---- measurement_analysis/ParseSphinxRefFWHMfile.m | 32 ++++++++++ measurement_analysis/Spots_LoadRef.m | 12 +++- 4 files changed, 112 insertions(+), 21 deletions(-) create mode 100644 measurement_analysis/ParseSphinxRefFWHMfile.m diff --git a/DisplaySpots.m b/DisplaySpots.m index 2925e3d..a872889 100644 --- a/DisplaySpots.m +++ b/DisplaySpots.m @@ -35,10 +35,9 @@ config=ConfigCheck(config,nSets,"config"); %% clear variables -[SP_Eks SP_Mms SP_FWHMs]=... - deal(missing(),missing(),missing()); -[MP_Eks MP_Mms MP_FWHMs MP_myPos]=... - deal(missing(),missing(),missing(),missing()); +[SP_Eks SP_Mms SP_FWHMs]=deal(missing(),missing(),missing()); +[MP_Eks MP_Mms MP_FWHMs MP_myPos]=deal(missing(),missing(),missing(),missing()); +[SPHINX_Mms SPHINX_FWHMs]=deal(missing(),missing()); %% get values @@ -63,8 +62,20 @@ MP_myPos=ExpandMat(MP_myPos,tmpMyPos'); end +% - acquire data from SPHINX +for iSet=1:nSets + % - parse DB file from SPHINX + [~,tmpSPHINX_Mms,tmpSPHINX_FWHMs]=Spots_LoadRef(machine(iSet),beamPart(iSet),config(iSet),"sphinx"); + % - store data + SPHINX_Mms=ExpandMat(SPHINX_Mms,repmat(tmpSPHINX_Mms,[1 size(tmpSPHINX_FWHMs,2)])); + SPHINX_FWHMs=ExpandMat(SPHINX_FWHMs,tmpSPHINX_FWHMs); +end + if (size(MP_FWHMs,3)~=size(SP_FWHMs,3)) - error("something wrong at parsing: %d MedPhys curves, and %d SteerPaz curves",size(MP_FWHMs,2),size(SP_FWHMs,3)); + error("something wrong at parsing: %d MedPhys curves, and %d SteerPaz curves",size(MP_FWHMs,3),size(SP_FWHMs,3)); +end +if (size(SPHINX_FWHMs,3)~=size(SP_FWHMs,3)) + error("something wrong at parsing: %d Sphinx curves, and %d SteerPaz curves",size(SPHINX_FWHMs,3),size(SP_FWHMs,3)); end %% show data @@ -84,7 +95,13 @@ xRef=MP_Mms(:,iISO,:); yRef=MP_FWHMs(:,iISO,:); yRef=Spots_FWHMRefSeries(yRef); ShowSeries(SP_Mms,SP_FWHMs,"R_{H_2O} [mm]","FWHM [mm]",["FWHM_x" "FWHM_y"],myTitles,yRef,xRef,["MP" "MP+" "MP-"]); -% - Medical Physics: various positions (grouped by line) +% - SPHINX vs Medical Physics +myTitles=strings(nDataSets,1); myTitles(:)=LineLabs; myTitles(end+1)="SPHINX vs MedPhys"; +iISO=(MP_myPos(:,1)==713E-3); +xRef=MP_Mms(:,iISO,:); yRef=MP_FWHMs(:,iISO,:); yRef=Spots_FWHMRefSeries(yRef); +ShowSeries(SPHINX_Mms,SPHINX_FWHMs,"R_{H_2O} [mm]","FWHM [mm]",["SPHINX_x","SPHINX_y","SPHINX_{mean}"],myTitles,yRef,xRef,["MP" "MP+" "MP-"]); + +%% - Medical Physics: various positions (grouped by line) myTitles=strings(nDataSets,1); myTitles(:)=LineLabs; myTitles(end)="MedPhys - various longitudinal positions"; myLeg=PosLabs; ShowSeries(MP_Mms,MP_FWHMs,"R_{H_2O} [mm]","FWHM [mm]",myLeg,myTitles); diff --git a/general/ShowSeries.m b/general/ShowSeries.m index c3edf64..410d018 100644 --- a/general/ShowSeries.m +++ b/general/ShowSeries.m @@ -1,4 +1,4 @@ -function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,refLeg,nCoupled) +function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,refLeg,pixY,pixX,pixLeg,nCoupled) % ShowStuff to display data series (xVals,yVals) as several plots in a single figure % % NOTA BENE: @@ -9,11 +9,17 @@ function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,r % - ref data can be: % . unique for all plots or a dedicated one for each plot -- size(refY,3); % . for each plot, a unique series (black) or a few (the first one black, the others grey) -- size(refY,2); +% - pix data can be: +% . unique for all plots or a dedicated one for each plot -- size(pixY,3); +% . for each plot, up to three series (squared magenta, pentagrammed cyan and hexagram green) -- size(pixY,2); if (~exist("myTitles","var")), myTitles=missing(); end if (~exist("refY","var")), refY=missing(); end if (~exist("refX","var")), refX=missing(); end if (~exist("refLeg","var")), refLeg=missing(); end + if (~exist("pixY","var")), pixY=missing(); end + if (~exist("pixX","var")), pixX=missing(); end + if (~exist("pixLeg","var")), pixLeg=missing(); end if (~exist("nCoupled","var")), nCoupled=1; end nSets=size(yVals,2); nShows=size(yVals,3); @@ -21,10 +27,6 @@ function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,r nSeriesLegendThresh=10; % max number of series before having a dedicated one for the legend figure(); - % - set grid of plots - nPlots=nShows; - if (all(~ismissing(myLeg)) && ( nPlots>nPlotsLegendThresh || nSets>nSeriesLegendThresh ) ), nPlots=nPlots+1; end - [nRows,nCols]=GetNrowsNcols(nPlots,nCoupled); % - set markers if ( nSets==2 ) markers=[ "o" "*" ]; @@ -38,6 +40,15 @@ function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,r refColors=NaN(size(refY,2),3); refColors(1,:)=[0 0 0]; if (size(refY,2)>1), refColors(2:end,:)=repmat([0.7 0.7 0.7],[size(refColors,1)-1 1]); end end + if (~ismissing(pixY)) + if ( size(pixY,2)<=3 ) + pixMarkers=[ "square" "pentagram" "hexagram" ]; + pixCm=[1 0 1 ; 0 1 1 ; 0 1 0]; % magenta, cyan and green + else + pixMarkers=strings(size(pixY,2),1); pixMarkers(:)="*"; + pixCm=colormap(jet(size(pixY,2))); + end + end % - legend tmpLeg=[]; if (all(~ismissing(refLeg))) @@ -46,11 +57,19 @@ function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,r if (all(~ismissing(myLeg))) tmpLeg=[tmpLeg myLeg]; end + if (all(~ismissing(pixLeg))) + tmpLeg=[tmpLeg pixLeg]; + end + nLeg=length(tmpLeg); + % - set grid of plots + nPlots=nShows; + if (~isempty(tmpLeg) && ( nShows>nPlotsLegendThresh || nLeg>nSeriesLegendThresh ) ), nPlots=nPlots+1; end + [nRows,nCols]=GetNrowsNcols(nPlots,nCoupled); % - actually plot for iPlot=1:nShows axs(iPlot)=subplot(nRows,nCols,iPlot); lFirst=true; - if (~ismissing(refX)) + if (~ismissing(refY)) for jRef=1:size(refY,2) if (lFirst), lFirst=false; else, hold on; end if (size(refX,2)==1), jX=1; else, jX=jRef; end @@ -67,21 +86,31 @@ function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,r if (size(xVals,3)==1), kX=1; else, kX=iPlot; end plot(xVals(:,jX,kX),yVals(:,iSet,iPlot),"-","Marker",markers(iSet),"Color",cm(iSet,:)); end + if (~ismissing(pixY)) + for jPix=1:size(pixY,2) + if (lFirst), lFirst=false; else, hold on; end + if (size(pixX,2)==1), jX=1; else, jX=jPix; end + if (size(pixX,3)==1), kX=1; else, kX=iPlot; end + if (size(pixY,2)==1), jY=1; else, jY=jPix; end + if (size(pixY,3)==1), kY=1; else, kY=iPlot; end + plot(pixX(:,jX,kX),pixY(:,jY,kY),"LineStyle","--","Marker",pixMarkers(jPix),"Color",pixCm(jPix,:)); + end + lFirst=false; + end if (length(showXlabel)==nShows), iX=iPlot; else, iX=1; end; xlabel(showXlabel(iX)); if (length(showYlabel)==nShows), iY=iPlot; else, iY=1; end; ylabel(showYlabel(iY)); grid on; if (all(~ismissing(myTitles))), title(myTitles(iPlot)); end - if (~isempty(tmpLeg) && ( nPlots<=nPlotsLegendThresh && nSets<=nSeriesLegendThresh ) ), legend(tmpLeg); end + if (~isempty(tmpLeg) && ( nPlots<=nPlotsLegendThresh && nLeg<=nSeriesLegendThresh ) ), legend(tmpLeg); end end % - legend plot - if (~isempty(tmpLeg) && ( nPlots>nPlotsLegendThresh || nSets>nSeriesLegendThresh ) ) + if (~isempty(tmpLeg) && ( nPlots>nPlotsLegendThresh || nLeg>nSeriesLegendThresh ) ) subplot(nRows,nCols,nPlots); lFirst=true; if (all(~ismissing(refLeg))) - plot(NaN(),NaN(),"k-"); - if (size(refY,2)>1) - hold on; plot(NaN(),NaN(),"-","Color",[.7 .7 .7]); - hold on; plot(NaN(),NaN(),"-","Color",[.7 .7 .7]); + for jRef=1:size(refY,2) + if (lFirst), lFirst=false; else, hold on; end + plot(NaN(),NaN(),refMarkers(jRef),"Color",refColors(jRef,:)); end lFirst=false; end @@ -89,6 +118,13 @@ function ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg,myTitles,refY,refX,r if (lFirst), lFirst=false; else, hold on; end plot(NaN(),NaN(),"-","Marker",markers(iSet),"Color",cm(iSet,:)); end + if (all(~ismissing(pixLeg))) + for jPix=1:size(pixY,2) + if (lFirst), lFirst=false; else, hold on; end + plot(NaN(),NaN(),"LineStyle","--","Marker",pixMarkers(jPix),"Color",pixCm(jPix,:)); + end + lFirst=false; + end legend(tmpLeg,"Location","best"); title("legend"); end % - general stuff diff --git a/measurement_analysis/ParseSphinxRefFWHMfile.m b/measurement_analysis/ParseSphinxRefFWHMfile.m new file mode 100644 index 0000000..c7b0311 --- /dev/null +++ b/measurement_analysis/ParseSphinxRefFWHMfile.m @@ -0,0 +1,32 @@ +function [Mms,FWHMs]=ParseSphinxRefFWHMfile(fileName,machine,beamPart) +% output vars: +% - size(FWHMs,1)=length(Mms); +% - size(FWHMs,2)=3 (ie FWHMx, FWHMy, mean(FWHMx,FWHMy)); + if (~exist("machine","var")), machine=missing(); end + if (~exist("beamPart","var")), beamPart=missing(); end + if (ismissing(machine)), machine="Sala1"; end % default + if (ismissing(beamPart)), beamPart="Prot"; end % default + fprintf("parsing Sphinx ref file for FWHM values %s ...\n",fileName); + switch upper(beamPart) + case {"P","PROT","PROTON"} + myData=readmatrix(fileName,"Range","A3:M11"); + case {"C","CARB","CARBON"} + myData=readmatrix(fileName,"Range","A15:M24"); + otherwise + error("Unknown particle %s!",beamPart); + end + Mms=myData(:,1); + switch upper(machine) + case "SALA1" + FWHMs=myData(:,2:4); + case "SALA2H" + FWHMs=myData(:,5:7); + case "SALA2V" + FWHMs=myData(:,11:13); + case "SALA3" + FWHMs=myData(:,8:10); + otherwise + error("Unknown machine: %s",machine); + end + fprintf("...data loaded: %d energies, %d series;\n",size(FWHMs,1),size(FWHMs,2)); +end diff --git a/measurement_analysis/Spots_LoadRef.m b/measurement_analysis/Spots_LoadRef.m index 34c013c..dd6c268 100644 --- a/measurement_analysis/Spots_LoadRef.m +++ b/measurement_analysis/Spots_LoadRef.m @@ -4,17 +4,23 @@ if (~exist("loc","var")), loc="ISO"; end myConfig=sprintf("%s,%s,%s",machine,beamPart,config); FullFileName=ReturnDefFile(spotType,myConfig); + [dMachine,dBeamPart,dFocus,dConfig]=DecodeConfig(myConfig); switch upper(spotType) case {"STEERINGPAZIENTI","SP","STEERPAZ"} [cyProgs,cyCodes,BARs,FWHMs,ASYMs,INTs]=ParseBeamProfileSummaryFiles(FullFileName,"CAM"); - if (length(cyProgs)<=1), error("...no data aquired!"); end + if (length(cyProgs)<1), error("...no data aquired!"); end % - convert cyCodes [Eks,~]=MapCyCodes(cyCodes,"Ek"); [Mms,~]=MapCyCodes(cyCodes,"Range"); pos=NaN(); case {"MEDICALPHYSICS","MP","MEDPHYS"} - [Eks,Mms,FWHMs,pos]=ParseMPrefFWHMfile(FullFileName,beamPart,loc); - if (length(Eks)<=1), error("...no data aquired!"); end + [Eks,Mms,FWHMs,pos]=ParseMPrefFWHMfile(FullFileName,dBeamPart,loc); + if (length(Eks)<1), error("...no data aquired!"); end + case "SPHINX" + [Mms,FWHMs]=ParseSphinxRefFWHMfile(FullFileName,dMachine,dBeamPart); + if (length(Mms)<1), error("...no data aquired!"); end + pos=NaN(); + Eks=NaN(); otherwise error("Unknown spot type %s!",spotType); end