dresp_dke_yp

       Calculation of the energy response function of detector (scintillator
       coupled with a photo-multiplier or diode). Compton scattering, photoelectric
       absorption and detector resolution due to counting noise are 
       taken into account.

       Input:

               - res: fit parameters of the energy resolution of the detectors [m,2]
               - phf: photofraction determined using the "mcdet.f" Monte-carlo hard x-ray absorption code at energies kd [m,p]
               - kd: reference energies for the photofraction [m,p]
               - k: photon energy (keV) [1,n]

       Output:
       
               - gkE: full detector's response matrix (keV-1) [n,n]


by Y.PEYSSON CEA-DRFC <yves.peysson@cea.fr>

0001 function gkE = dresp_dke_yp(res,phf,kd,k)
0002 %
0003 %       Calculation of the energy response function of detector (scintillator
0004 %       coupled with a photo-multiplier or diode). Compton scattering, photoelectric
0005 %       absorption and detector resolution due to counting noise are
0006 %       taken into account.
0007 %
0008 %       Input:
0009 %
0010 %               - res: fit parameters of the energy resolution of the detectors [m,2]
0011 %               - phf: photofraction determined using the "mcdet.f" Monte-carlo hard x-ray absorption code at energies kd [m,p]
0012 %               - kd: reference energies for the photofraction [m,p]
0013 %               - k: photon energy (keV) [1,n]
0014 %
0015 %       Output:
0016 %
0017 %               - gkE: full detector's response matrix (keV-1) [n,n]
0018 %
0019 %
0020 %by Y.PEYSSON CEA-DRFC <yves.peysson@cea.fr>
0021 %
0022 if nargin < 4,
0023    error('Wrong number of input arguments for dresp_dke_yp');
0024 end
0025 %
0026 [qe,me,mp,mn,e0,mu0,re,mc2,clum,alpha] = pc_dke_yp;%Physics constant
0027 %
0028 if min(size(k)) > 1,
0029         error('k must be a vector');
0030 end
0031 k = k(:);
0032 k = k';
0033 nk = length(k);
0034 deltak = k(2)-k(1);
0035 E = k;
0036 %
0037 gkE = zeros(nk,nk);
0038 %
0039 Rk = (interp1(kd,phf,k,'spline'))';
0040 dEk = (res(1)*k.^res(2)).*k;
0041 Ec = (2*k.^2)./(2*k+mc2);
0042 dEc = (res(1)*Ec.^res(2)).*Ec;
0043 %
0044 for i=1:nk,
0045         ic = max(find(E<=Ec(i)));
0046         if isempty(ic),ic = 0;end
0047         %
0048         if dEk(i)~=0,
0049                 gkEph = Rk(i)*2*sqrt(log(2)/pi)*exp(-4*log(2)*(E - k(i)).^2/dEk(i)^2)/dEk(i);
0050                 if ic==0,
0051                         gkE(i,:) = gkEph;%/Rk(i);
0052                 else
0053                         E1 = E(1:ic);
0054                         E2 = E(ic+1:nk);
0055                         gkEco1 = (1 - Rk(i))/(0.5*dEc(i)/2/sqrt(log(2)/pi) + Ec(i))*ones(1,size(E1,2));
0056                         gkEco2 = (1 - Rk(i))/(0.5*dEc(i)/2/sqrt(log(2)/pi) + Ec(i))*exp(-4*log(2)*(E2 - E2(1)).^2/dEc(i)^2);
0057                         gkEco = [gkEco1,gkEco2];
0058                         gkE(i,:) = gkEph + gkEco;
0059                 end
0060         else
0061                 gkEph0 = speye(i,nk);
0062                 gkEph = Rk(i)*gkEph0/deltak;
0063                 if ic==0,
0064                         gkE(i,:) = gkEph;%/Rk(i);
0065                 else
0066                         E1 = E(1:ic);
0067                         E2 = E(ic+1:nk);                    
0068                         gkEco1 = ((1 - Rk(i))/Ec(i))*ones(1,size(E1,2));
0069                         gkEco2 = 0.0*ones(1,size(E2,2));
0070                         gkEco = [gkEco1,gkEco2];
0071                         gkE(i,:) = gkEph + gkEco;
0072                 end
0073         end
0074 end
0075 
0076