diff --git a/CreateQAfiles.m b/CreateQAfiles.m index d63d59a..6f27e49 100644 --- a/CreateQAfiles.m +++ b/CreateQAfiles.m @@ -14,27 +14,28 @@ pathToLibrary=".\"; addpath(genpath(pathToLibrary)); end +% - include MachineRefs +pathToLibrary="../MachineRefs"; +addpath(genpath(pathToLibrary)); %% settings % ------------------------------------------------------------------------- % USER's input data -kPath="S:\Accelerating-System\Accelerator-data"; -% kPath="K:"; - -beamPart="CARBON"; machine="ISO3"; +beamPart="PROTON"; 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); @@ -48,20 +49,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"; % ------------------------------------------------------------------------- @@ -72,8 +75,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)) @@ -88,7 +91,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)); 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/DisplayLGENs.m b/DisplayLGENs.m new file mode 100644 index 0000000..afd82cc --- /dev/null +++ b/DisplayLGENs.m @@ -0,0 +1,81 @@ +% {}~ + +%% 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=["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=["PROTON" "CARBON" ]; +% machine="ISO2"; +% beamPart=["PROTON" "CARBON" ]; +config="TM"; % select configuration: TM, RFKO +myTitle="LGEN values"; +% ------------------------------------------------------------------------- + +% ------------------------------------------------------------------------- +% 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 +[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,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,"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 diff --git a/DisplayRefValues.m b/DisplayRefValues.m new file mode 100644 index 0000000..5f4fa99 --- /dev/null +++ b/DisplayRefValues.m @@ -0,0 +1,63 @@ +% {}~ + +%% 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=ConfigCheck(beamPart,nSets,"beamPart"); +machine=ConfigCheck(machine,nSets,"machine"); +mSets=max([length(showX) length(showXlabel) length(showY) length(showYlabel)]); +showX=ConfigCheck(showX,mSets,"showX"); +showXlabel=ConfigCheck(showXlabel,mSets,"showXlabel"); +showY=ConfigCheck(showY,mSets,"showY"); +showYlabel=ConfigCheck(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]=ExtractFromTable(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 +ShowSeries(xVals,yVals,showXlabel,showYlabel,myLeg); + +%% local functions diff --git a/DisplaySpots.m b/DisplaySpots.m new file mode 100644 index 0000000..a872889 --- /dev/null +++ b/DisplaySpots.m @@ -0,0 +1,119 @@ +% {}~ + +%% description +% this is a script which displays beam spots; + +%% include libraries +% - include Matlab libraries +if (~exist("pathToLibrary","var")) + pathToLibrary="./"; + addpath(genpath(pathToLibrary)); + pathToLibrary="../MachineRefs"; + addpath(genpath(pathToLibrary)); +end + +%% clear +clear all; +close all; + +%% settings + +% ------------------------------------------------------------------------- +% USER's input data +machine=["Sala1" "Sala2H" "Sala2V" "Sala3"]; +beamPart="PROTON"; +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()); +[SPHINX_Mms SPHINX_FWHMs]=deal(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 + +% - 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,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 +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-"]); + +% - 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); + +% - 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/MeVvsCyCo_P.xlsx b/MeVvsCyCo_P.xlsx deleted file mode 100644 index 7fa3748..0000000 Binary files a/MeVvsCyCo_P.xlsx and /dev/null differ diff --git a/PSmapping.xlsx b/PSmapping.xlsx deleted file mode 100644 index 4df5af4..0000000 Binary files a/PSmapping.xlsx and /dev/null differ diff --git a/README.md b/README.md index 7822937..1a05841 100644 --- a/README.md +++ b/README.md @@ -8,6 +8,7 @@ For the time being: | folder name | description | |-----|-----| +| `displayLib` | library of functions used by `Display*` main scripts | | `educational` | scripts of general interest showing quantities relevant to CNAO (e.g. bar of charge, beam quantities vs energy, etc...) | | `general` | scripts of general interest (e.g. math, MatLab environment, etc...) | | `MADX-optics` | scripts for plotting optics computed by MADX | diff --git a/displayLib/ConfigCheck.m b/displayLib/ConfigCheck.m new file mode 100644 index 0000000..bddd5e9 --- /dev/null +++ b/displayLib/ConfigCheck.m @@ -0,0 +1,30 @@ +function OutVar=ConfigCheck(InVar,nSets,myName) +% ConfigCheck(InVar,nSets,myName) consistency check of user input data; +% the function simply checks that the +% array InVar has nSets non-missing +% values; +% in case nSets>1 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 (nIn1 && 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 401c88d..0f8814f 100644 --- a/educational/ShowMCS.m +++ b/educational/ShowMCS.m @@ -13,51 +13,49 @@ 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" + +% 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=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 -run(".\particleData.m"); +%% Load particle data +% returns: myM [MeV/c2], myEk [MeV], myZ [], unitEk ("MeV" for protons, "MeV/u" for others); +run(".\setParticle.m"); +%% 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)) +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_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,unitEk)); set(gca, 'YScale', 'log'); end diff --git a/educational/particleData.m b/educational/particleData.m index be2be55..d3f3585 100644 --- a/educational/particleData.m +++ b/educational/particleData.m @@ -2,12 +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 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..410d018 --- /dev/null +++ b/general/ShowSeries.m @@ -0,0 +1,135 @@ +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: +% - 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); +% - 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); + 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 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 + 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))) + tmpLeg=refLeg; + end + 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(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 + 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 (~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 && nLeg<=nSeriesLegendThresh ) ), legend(tmpLeg); end + end + % - legend plot + if (~isempty(tmpLeg) && ( nPlots>nPlotsLegendThresh || nLeg>nSeriesLegendThresh ) ) + subplot(nRows,nCols,nPlots); + lFirst=true; + if (all(~ismissing(refLeg))) + 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 + for iSet=1:nSets + 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 + if (all(~ismissing(myTitles))) + if (length(myTitles)>nShows), sgtitle(myTitles(end)); end + end + % linkaxes(axs,"xy"); +end 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; 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/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_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..f181a53 --- /dev/null +++ b/measurement_analysis/Spots_LoadEffectiveSpecs.m @@ -0,0 +1,36 @@ +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 + [tmpEks,tmpMms,tmpFWHMs]=Spots_LoadRef(machine(iSet),beamPart,config,spotType); + % - store data + FWHMs=ExpandMat(FWHMs,tmpFWHMs); + 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 + 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 diff --git a/measurement_analysis/Spots_LoadRef.m b/measurement_analysis/Spots_LoadRef.m new file mode 100644 index 0000000..dd6c268 --- /dev/null +++ b/measurement_analysis/Spots_LoadRef.m @@ -0,0 +1,27 @@ +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); + [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 + % - convert cyCodes + [Eks,~]=MapCyCodes(cyCodes,"Ek"); + [Mms,~]=MapCyCodes(cyCodes,"Range"); + pos=NaN(); + case {"MEDICALPHYSICS","MP","MEDPHYS"} + [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 +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/AcquireLGENValues.m b/operations/AcquireLGENValues.m index 1df8d2b..b6efde2 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,159 +22,11 @@ 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 - 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)) ; % - 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 + %% do the job + 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/DecodeConfig.m b/operations/DecodeConfig.m new file mode 100644 index 0000000..9cebff2 --- /dev/null +++ b/operations/DecodeConfig.m @@ -0,0 +1,101 @@ +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"; focus=""; + % - actual requests by user + if (~ismissing(myConfig) && strlength(myConfig)>0) + myInfo=split(myConfig,","); + if (length(myInfo)>=1) + if ( strlength(myInfo(1))>0 ), machine=myInfo(1); end + end + if (length(myInfo)>=2) + if ( strlength(myInfo(2))>0 ), beamPart=myInfo(2); end + end + if (length(myInfo)>=3) + if ( strlength(myInfo(3))>0 ), config=myInfo(3); end + end + if (length(myInfo)>=4) + 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"} + 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") || strlength(focus)==0), focus="FG"; end + case {"CARBON","CARBONIO","CARB","C"} + beamPart="C"; + if (~exist("focus","var") || strlength(focus)==0), 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 + if (strlength(focus)>0) + error("focus %s NOT recognised!",focus); + end + end + if (strcmp(machine,"SYNCHRO")), focus=""; end + +end diff --git a/operations/LGENvisualCheck.m b/operations/LGENvisualCheck.m index ff635d3..93359a7 100644 --- a/operations/LGENvisualCheck.m +++ b/operations/LGENvisualCheck.m @@ -1,23 +1,65 @@ -function LGENvisualCheck(psNames,xVals,xLab,yVals,yLab,magNames,myTitle,LGENs) - if (~exist("LGENs","var")), LGENs=psNames; end +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("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(); - nSets=length(LGENs); - % automatically set grid of subplots - [nRows,nCols]=GetNrowsNcols(nSets); - - for iSet=1:nSets - subplot(nRows,nCols,iSet); - jj=find(strcmp(psNames,LGENs(iSet))); - if ( isempty(jj) ) - warning("...unable to find LGEN %s in list of PS names!",LGENs(iSet)); - continue + cm=colormap(turbo(nSets)); + % - set grid of plots + nPlots=nMagnets; + if (~ismissing(myLegend)), nPlots=nPlots+1; end + [nRows,nCols]=GetNrowsNcols(nPlots); + % - set markers + if ( nSets==2 ) + markers=[ "o" "*" ]; + else + markers=strings(nSets,1); + markers(:)="."; + end + % - actually plot + for iMagnet=1:nMagnets + subplot(nRows,nCols,iMagnet); + lFirst=true; + for iSet=1:nSets + jj=find(strcmp(magNames(:,iSet),magNames2show(iMagnet))); + if ( isempty(jj) ) + warning("...unable to find LGEN %s in list of PS names!",magNames2show(iMagnet)); + continue + end + if (lFirst), lFirst=false; else, hold on; end + [~,IDs]=sort(xVals(:,iSet)); + plot(xVals(IDs,iSet),yVals(IDs,jj,iSet),"-","Marker",markers(iSet),"Color",cm(iSet,:)); end - plot(xVals,yVals(:,jj),'*'); - xlabel(xLab); ylabel(yLab); - grid on; - title(sprintf("%s (aka %s)",magNames(jj),LGENs(iSet))); + xlabel(xLab); ylabel(yLab); grid on; + title(sprintf("%s (aka %s)",magNames2show(iMagnet),psNames(jj,iSet))); end - - sgtitle(myTitle); + % - legend plot + if (~ismissing(myLegend)) + 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(myLegend); + end + % - general stuff + if (~ismissing(myTitle)), sgtitle(myTitle); end end 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 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 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