R5-X2 - Calculation of the bremsstrahlung emission as measured by the HXR diagnostic taking into account of the detector's response function.
Calculates the bremsstrahlung emission as measured by the HXR diagnostic
taking into account of the detector's response function.
Call first the function bremchord_dke_yp.m for calculating "vschords" used in line integration
Then call bremsstrahlung_dke_yp.m for calculating "brem_bd_out" used in the convolution
INPUTS:
- Zbremplasma.kphot: photon energies (keV) [1,nk] (given by bremsstrahlung_dke_yp.m)
- Zbremplasma.brem_bd_out: Plasma photon number per unit energy (keV) and time (s) for each chord [nchord,nk] (given by bremsstrahlung_dke_yp.m)
- hxr: Bremsstrahlung structure
- hxrparam : Bremsstrahlung calculation and display parameter structure
- dt: integration time step (s) [1,1]
OUTPUTS:
- Zbremdiag: Diagnostic convolved bremstrahlung spectrum data structure
* Zbremdiag.kinterp: Photon energy spectrum grid for response function calculations (kphotmin:kphotmax) [1,1000]
* Zbremdiag.gkEinterp_out: Detector x-ray response function {nchord}
* Zbremdiag.kfit: Discrete pulse energy spectrum (keV) [1,nkfit]
* Zbremdiag.dkfit: Discrete pulse energy spectrum step width (keV) [1,nkfit]
* Zbremdiag.brem_bd_plasma: Plasma photon number per unit energy (keV) and time (s) for each chord [nchord,nkfit]
* Zbremdiag.brem_bd_diag: Pulse number per unit energy (keV) and time (s) for each chord [nchord,nkfit]
* Zbremdiag.tphot_plasma: Plasma photon temperature (keV) [1,nchord]
* Zbremdiag.afit_plasma: A(k=0) [s^-1]','k.dN_phot/dk.dt = A.exp(-k/T_phot)@k=0 [50-110 keV] [1,nchord]
* Zbremdiag.tphot_plasma_exp: Pulse photon temperature (keV) [1,nchord]
* Zbremdiag.afit_plasma_exp: A*(E=0) [s-1]','E.dN_pulse/dE.dt.eta_E = A.exp(-E/T_phot*)@E=0 [50-110 keV] [1,nchord]
* Zbremdiag.meanabsx: Mean stopping efficiency [1,nchord]
by Joan Decker (CEA/IRFM,joan.decker@cea.fr) and Yves Peysson(CEA/IRFM,yves.peysson@cea.fr)0001 function [Zbremdiag] = bremdiag_dke_yp(Zbremplasma,hxr,hxrparam,dt) 0002 % 0003 % R5-X2 - Calculation of the bremsstrahlung emission as measured by the HXR diagnostic taking into account of the detector's response function. 0004 % 0005 %Calculates the bremsstrahlung emission as measured by the HXR diagnostic 0006 %taking into account of the detector's response function. 0007 %Call first the function bremchord_dke_yp.m for calculating "vschords" used in line integration 0008 %Then call bremsstrahlung_dke_yp.m for calculating "brem_bd_out" used in the convolution 0009 % 0010 % INPUTS: 0011 % 0012 % - Zbremplasma.kphot: photon energies (keV) [1,nk] (given by bremsstrahlung_dke_yp.m) 0013 % - Zbremplasma.brem_bd_out: Plasma photon number per unit energy (keV) and time (s) for each chord [nchord,nk] (given by bremsstrahlung_dke_yp.m) 0014 % - hxr: Bremsstrahlung structure 0015 % - hxrparam : Bremsstrahlung calculation and display parameter structure 0016 % - dt: integration time step (s) [1,1] 0017 % 0018 % OUTPUTS: 0019 % 0020 % - Zbremdiag: Diagnostic convolved bremstrahlung spectrum data structure 0021 % * Zbremdiag.kinterp: Photon energy spectrum grid for response function calculations (kphotmin:kphotmax) [1,1000] 0022 % * Zbremdiag.gkEinterp_out: Detector x-ray response function {nchord} 0023 % * Zbremdiag.kfit: Discrete pulse energy spectrum (keV) [1,nkfit] 0024 % * Zbremdiag.dkfit: Discrete pulse energy spectrum step width (keV) [1,nkfit] 0025 % * Zbremdiag.brem_bd_plasma: Plasma photon number per unit energy (keV) and time (s) for each chord [nchord,nkfit] 0026 % * Zbremdiag.brem_bd_diag: Pulse number per unit energy (keV) and time (s) for each chord [nchord,nkfit] 0027 % * Zbremdiag.tphot_plasma: Plasma photon temperature (keV) [1,nchord] 0028 % * Zbremdiag.afit_plasma: A(k=0) [s^-1]','k.dN_phot/dk.dt = A.exp(-k/T_phot)@k=0 [50-110 keV] [1,nchord] 0029 % * Zbremdiag.tphot_plasma_exp: Pulse photon temperature (keV) [1,nchord] 0030 % * Zbremdiag.afit_plasma_exp: A*(E=0) [s-1]','E.dN_pulse/dE.dt.eta_E = A.exp(-E/T_phot*)@E=0 [50-110 keV] [1,nchord] 0031 % * Zbremdiag.meanabsx: Mean stopping efficiency [1,nchord] 0032 % 0033 % by Joan Decker (CEA/IRFM,joan.decker@cea.fr) and Yves Peysson(CEA/IRFM,yves.peysson@cea.fr) 0034 % 0035 if nargin < 4, 0036 error('Not enough arguments in bremshtrahlung_dke_yp') 0037 end 0038 % 0039 method = 'cubic'; 0040 % 0041 kphot = Zbremplasma.kphot; 0042 brem_bd_in = Zbremplasma.brem_bd_out; 0043 % 0044 kphotmin = min(kphot);%hxr.kdiag_hxr(1,1); 0045 kphotmax = max(kphot);%hxr.kdiag_hxr(2,end-1); 0046 % 0047 % kfit = (hxr.kdiag_hxr(1,1:end-1) + hxr.kdiag_hxr(2,1:end-1))/2;%Measured photon energy spectrum (keV) [1,k] 0048 % kfit_S = [3*kphot(1) - kphot(2),kphot(1:end-1) + kphot(2:end),3*kphot(end) - kphot(end-1)]/2; 0049 % kfit = kphot;%Measured photon energy spectrum (keV) [1,k] 0050 % 0051 % kfit = hxr.kfit;%Measured photon energy spectrum (keV) [1,k] 0052 % kfit_S = [3*kfit(1) - kfit(2),kfit(1:end-1) + kfit(2:end),3*kfit(end) - kfit(end-1)]/2; 0053 % 0054 kfit_S = hxr.kfit;%Measured photon energy spectrum (keV) [1,k] 0055 kfit = (kfit_S(1:end-1) + kfit_S(2:end))/2; 0056 % 0057 % ikfitmin = find(kfit >= hxr.kdiag_hxr(1,end),1,'first'); 0058 % ikfitmax = find(kfit <= hxr.kdiag_hxr(2,end),1,'last'); 0059 ikfitmin = find(kfit >= hxrparam.kmin,1,'first'); 0060 ikfitmax = find(kfit <= hxrparam.kmax,1,'last'); 0061 %kfitmin = kfit(ikfitmin); 0062 %kfitmax = kfit(ikfitmax); 0063 % dkfit = hxr.kdiag_hxr(2,1:end-1) - hxr.kdiag_hxr(1,1:end-1);%Channel width of the measured photon energy spectrum (keV) [1,k] 0064 dkfit = diff(kfit_S);%Channel width of the measured photon energy spectrum (keV) [1,k] 0065 nkfit = length(dkfit); 0066 % 0067 kinterp = linspace(kphotmin,kphotmax,1001);%Uniform grid only 0068 dkinterp = kinterp(2) - kinterp(1); 0069 brem_bd_interp = interp1(kphot,brem_bd_in.',kinterp,method).'; 0070 % 0071 if size(brem_bd_in,1) == 1, 0072 brem_bd_interp = brem_bd_interp.'; 0073 end 0074 % 0075 for ik = 1:nkfit, 0076 brem_bd_plasma(:,ik) = sum(brem_bd_interp(:,find(kinterp >= kfit_S(ik),1,'first'):find(kinterp < kfit_S(ik+1),1,'last')).').'*dkinterp;%number of photons in energy channel ik per unit time 0077 end 0078 brem_bd_plasma = abs(brem_bd_plasma); 0079 % 0080 for ic = 1:size(brem_bd_plasma,1), 0081 [Tph,eTph,A,eA,Q] = tphfit_dke_yp(kfit(ikfitmin:ikfitmax)',dt*brem_bd_plasma(ic,ikfitmin:ikfitmax)',sqrt(dt*brem_bd_plasma(ic,ikfitmin:ikfitmax))'); 0082 tphot_plasma(ic) = Tph; 0083 etphot_plasma(ic) = eTph; 0084 afit_plasma(ic) = A; 0085 eafit_plasma(ic) = eA; 0086 Q_plasma(ic) = Q; 0087 end 0088 % 0089 %wbrem_bd_plasma = brem_bd_plasma.*repmat(kfit,[size(brem_bd_plasma,1),1])./repmat(dkfit,[size(brem_bd_plasma,1),1]);% per keV but taking into account of the channel energy width 0090 %for ic = 1:size(brem_bd_plasma,1), 0091 % pfit = polyfit(kfit(ikfitmin:ikfitmax),log(wbrem_bd_plasma(ic,ikfitmin:ikfitmax)),1); 0092 % tphot_plasma(ic) = -1/pfit(1); 0093 % afit_plasma(ic) = exp(pfit(2)); 0094 %end 0095 % 0096 % Convolution with detector response 0097 % 0098 kd = hxr.kd_hxr;%reference energies for the photofraction [m,p] 0099 phf = hxr.phf_hxr;%photofraction determined using the "MCDET.f" Monte-carlo hard x-ray absorption code at energies kd [m,p] 0100 res = hxr.res_hxr;%fit parameters of the energy resolution of the detectors [m,2] 0101 ef = hxr.ef_hxr;%thickness of the vacuum window (mm) [1,m] 0102 tf = hxr.tf_hxr;%type of material for the vacuum window (Al ou Be) [1,m] 0103 ea = hxr.ea_hxr;%thickness of the absorber (mm) [1,m] 0104 ta = hxr.ta_hxr;%type of material for the absorber (Al,Fe,Ge,Be ou Pb) [1,m] 0105 ec = hxr.ec_hxr;%thickness of the detector shield (mm) [1,m] 0106 tc = hxr.tc_hxr;%type of material for the detector shield (Al ou Be) [1,m] 0107 ed = hxr.ed_hxr;%thickness of the detector (mm) [1,m] 0108 td = hxr.td_hxr;%type of material for the detector (BGO,CsI,NaI,Ge or CdTe) [1,m] 0109 % 0110 meanabsx = zeros(size(brem_bd_plasma));%Mean attenuation per channel 0111 brem_bd_diag_interp = zeros(size(brem_bd_interp)); 0112 brem_bd_diag = zeros(size(brem_bd_plasma)); 0113 % 0114 for ic = 1:size(brem_bd_interp,1), 0115 if ic == 1, 0116 gkEinterp1 = dresp_dke_yp(res(ic,:),phf(ic,:),kd(ic,:),kinterp); 0117 normgkEinterp1 = trapz(kinterp,gkEinterp1');%To enforce that the norm is unity 0118 gkEinterp = gkEinterp1;%(gkEinterp1'./(normgkEinterp1(:)*ones(1,length(kinterp)))')'; 0119 % plot(kinterp,trapz(kinterp,gkEinterp1') + gkEinterp1(:,1)'*kinterp(1));%test of the norm 0120 % plot(kinterp,trapz(kinterp,gkEinterp'));%test of the norm 0121 % 0122 absxinterp = absx_dke_yp(kinterp,0,0,ef(ic),tf(ic,:),ea(ic),ta(ic,:),ec(ic),tc(ic,:),ed(ic),td(ic,:));%Air contribution is neglected (k > 10 keV only) 0123 Tabsxinterp = repmat(absxinterp(:),[1,size(brem_bd_interp,2)]); 0124 else 0125 if (kd(ic-1) ~= kd(ic)) | (res(ic-1) ~= res(ic)) | (phf(ic-1) ~= phf(ic)) 0126 gkEinterp1 = dresp_dke_yp(res(ic,:),phf(ic,:),kd(ic,:),ic,kinterp);%Recalculated only when detector type has changed concerning the response function 0127 normgkEinterp1 = trapz(kinterp,gkEinterp1');%To enforce that the norm is unity 0128 gkEinterp = (gkEinterp1'./(normgkEinterp1(:)*ones(1,length(kinterp)))')'; 0129 end 0130 % 0131 if (ef(ic-1) ~= ef(ic)) | (tf(ic-1) ~= tf(ic)) | (ea(ic-1) ~= ea(ic)) | (ta(ic-1) ~= ta(ic)) | (ec(ic-1) ~= ec(ic)) | (tc(ic-1) ~= tc(ic)) | (ed(ic-1) ~= ed(ic)) | (td(ic-1) ~= td(ic)) 0132 absxinterp = absx_dke_yp(kinterp,0,0,ef(ic),tf(ic,:),ea(ic),ta(ic,:),ec(ic),tc(ic,:),ed(ic),td(ic,:));%Air contribution is neglected (k > 10 keV only) 0133 Tabsxinterp = repmat(absxinterp(:),[1,size(brem_bd_interp,2)]); 0134 end 0135 end 0136 % 0137 gkEinterp_out{ic} = gkEinterp;%For display of the response function 0138 % 0139 brem_bd_diag_interp(ic,:) = (((gkEinterp'.*Tabsxinterp')*dkinterp*(brem_bd_interp(ic,:))'))'; 0140 % brem_bd_diag_interp(ic,:) = (fix(((eye(size(gkEinterp')).*Tabsxinterp)*(brem_bd_interp(ic,:))')))';%detector response neglected 0141 % 0142 for ik = 1:nkfit, 0143 meanabsx(ic,ik) = sum(absxinterp(min(find(kinterp >= kfit_S(ik))):max(find(kinterp < kfit_S(ik+1))))')'/length(absxinterp(min(find(kinterp >= kfit_S(ik))):max(find(kinterp < kfit_S(ik+1)))).'); 0144 end 0145 % 0146 for ik = 1:nkfit, 0147 brem_bd_diag(ic,ik) = sum(brem_bd_diag_interp(ic,min(find(kinterp >= kfit_S(ik))):max(find(kinterp < kfit_S(ik+1)))).').'*dkinterp;%number of pulses in energy channel ik per unit time 0148 end 0149 % 0150 brem_bd_diag = abs(brem_bd_diag); 0151 % 0152 [Tph,eTph,A,eA,Q] = tphfit_dke_yp(kfit(ikfitmin:ikfitmax)',dt*brem_bd_diag(ic,ikfitmin:ikfitmax)',sqrt(dt*brem_bd_diag(ic,ikfitmin:ikfitmax))'); 0153 if Tph <=0, %Remove negative or null values 0154 Tph = NaN; 0155 eTph = NaN; 0156 A = NaN; 0157 eA = NaN; 0158 Q = NaN; 0159 end 0160 % 0161 tphot_plasma_exp(ic) = Tph; 0162 etphot_plasma_exp(ic) = eTph; 0163 afit_plasma_exp(ic) = A; 0164 eafit_plasma_exp(ic) = eA; 0165 Q_plasma_exp(ic) = Q; 0166 0167 % wbrem_bd_plasma_exp(ic,:) = brem_bd_diag(ic,:).*kfit./dkfit./meanabsx(ic,:);% per keV but taking into account of the channel energy width and mean attenuation per channel 0168 % pfit = polyfit(kfit(ikfitmin:ikfitmax),log(wbrem_bd_plasma_exp(ic,ikfitmin:ikfitmax)),1); 0169 % tphot_plasma_exp(ic) = -1/pfit(1); 0170 % afit_plasma_exp(ic) = exp(pfit(2)); 0171 % if tphot_plasma_exp(ic) <=0, %Remove negative or null values 0172 % tphot_plasma_exp(ic) = NaN; 0173 % afit_plasma_exp(ic) = NaN; 0174 % end 0175 % 0176 info_dke_yp(2,['Convolution with detector''s response function for chord #',int2str(ic),'/',int2str(size(brem_bd_interp,1)),' done']); 0177 end 0178 % 0179 Zbremdiag.kinterp = kinterp; 0180 Zbremdiag.gkEinterp_out = gkEinterp_out; 0181 Zbremdiag.kfit = kfit; 0182 Zbremdiag.dkfit = dkfit; 0183 Zbremdiag.brem_bd_plasma = brem_bd_plasma; 0184 Zbremdiag.brem_bd_diag = brem_bd_diag; 0185 Zbremdiag.tphot_plasma = tphot_plasma; 0186 Zbremdiag.etphot_plasma = etphot_plasma; 0187 Zbremdiag.afit_plasma = afit_plasma; 0188 Zbremdiag.eafit_plasma = eafit_plasma; 0189 Zbremdiag.Q_plasma = Q_plasma; 0190 Zbremdiag.tphot_plasma_exp = tphot_plasma_exp; 0191 Zbremdiag.etphot_plasma_exp = etphot_plasma_exp; 0192 Zbremdiag.afit_plasma_exp = afit_plasma_exp; 0193 Zbremdiag.eafit_plasma_exp = eafit_plasma_exp; 0194 Zbremdiag.Q_plasma_exp = Q_plasma_exp; 0195 Zbremdiag.meanabsx = meanabsx; 0196 Zbremdiag.dt = dt;