fpilup_dke_yp

    Study of of pile-up and dead-time for x-ray spectrometry at high counting rate
    Global study based on Poisson statistics. Monte-Carlo simulation.

    Input: 

        - diag: diagnostic structure (string) [1,m]
       - nreset: random number generator reset number
        - Ni0: number of pulses [1,1]
          (default = 10000)
        - kmin0: lower limit of the photon energy spectrum (keV) [1,1]
          (default = kmin (see the listing)
        - k0: photon reference energy (keV) [1,1]
          (default = 70)
        - ech: scaling factor to reduce Monte-Carlo calculation time [1,1]
          (default = 10)
        - typespec: (1) line spectrum, (2) exponential spectrum [1,1]
          (default = 2)
        - pur: (1) pile-up rejector on, or (0) off [1,1]
          (default = 1)

    Output: 

        - ki:
        - ko:
        - kADC:
        - kkADC:
        - emp:
        - kcode:
        - khisto: energy range (keV) [1,nenergy]
        - khisto1: energy range (keV) [1,nenergy]
        - histoki: input spectrum (plasma one) [1,nenergy]
        - histoko: output spectrum (without ADC dead time) [1,nenergy]
        - histokADC: output spectrum (with ADC dead time) [1,nenergy]
        - histokcode: output spectrum (with or without pur) [1,nenergy]


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

    Diagnostics parameters

0001 function [ki,ko,kADC,kkADC,emp,kcode,khisto,khisto1,histoki,histoko,histokADC,histokcode] = fpilupyp(diag,nreset,Ni0,k0min,k0,ech,typspec,pur)
0002 %
0003 %    Study of of pile-up and dead-time for x-ray spectrometry at high counting rate
0004 %    Global study based on Poisson statistics. Monte-Carlo simulation.
0005 %
0006 %    Input:
0007 %
0008 %        - diag: diagnostic structure (string) [1,m]
0009 %       - nreset: random number generator reset number
0010 %        - Ni0: number of pulses [1,1]
0011 %          (default = 10000)
0012 %        - kmin0: lower limit of the photon energy spectrum (keV) [1,1]
0013 %          (default = kmin (see the listing)
0014 %        - k0: photon reference energy (keV) [1,1]
0015 %          (default = 70)
0016 %        - ech: scaling factor to reduce Monte-Carlo calculation time [1,1]
0017 %          (default = 10)
0018 %        - typespec: (1) line spectrum, (2) exponential spectrum [1,1]
0019 %          (default = 2)
0020 %        - pur: (1) pile-up rejector on, or (0) off [1,1]
0021 %          (default = 1)
0022 %
0023 %    Output:
0024 %
0025 %        - ki:
0026 %        - ko:
0027 %        - kADC:
0028 %        - kkADC:
0029 %        - emp:
0030 %        - kcode:
0031 %        - khisto: energy range (keV) [1,nenergy]
0032 %        - khisto1: energy range (keV) [1,nenergy]
0033 %        - histoki: input spectrum (plasma one) [1,nenergy]
0034 %        - histoko: output spectrum (without ADC dead time) [1,nenergy]
0035 %        - histokADC: output spectrum (with ADC dead time) [1,nenergy]
0036 %        - histokcode: output spectrum (with or without pur) [1,nenergy]
0037 %
0038 %
0039 %by Y.PEYSSON CEA-DRFC <yves.peysson@cea.fr>
0040 %
0041 %    Diagnostics parameters
0042 %
0043 tok = diag.tok;
0044 namediag = diag.namediag;
0045 theta = diag.theta;
0046 tc0 = diag.tc0;
0047 tpeak = diag.tpeak; 
0048 tl0 = diag.tl0;
0049 tp0 = diag.tp0;
0050 beta = diag.beta;
0051 kmin = diag.kmin;
0052 kmin0 = diag.kmin0;
0053 kmax = diag.kmax;
0054 nhisto = diag.nhisto;
0055 %
0056 tp = tp0*theta*ech;
0057 tl = tl0*theta*ech;
0058 tc = tc0*ech;
0059 tpeak = tpeak*ech;
0060 Ni = Ni0/ech;
0061 %
0062 alpha = beta*tp/(tl-tp);
0063 %
0064 %    Energy spectrum
0065 %
0066 khisto = linspace(kmin,kmax,nhisto+1);
0067 dk = khisto(2) - khisto(1);
0068 khisto1 = linspace(kmin+dk/2,kmax-dk/2,nhisto);
0069 %
0070 if typspec==1,
0071      ki = k0*ones(Ni,1);
0072     [histoki,dummy] = hist_dke_yp(ki,khisto);
0073     krefi = ki;
0074 elseif typspec==2,
0075      ki = kmin0-k0*log(rand(Ni,1));
0076      [histoki,dummy] = hist_dke_yp(ki,khisto);
0077      pfit = polyfit(khisto1,log(histoki),1);
0078      krefi = -1/pfit(1);afiti = exp(pfit(2));
0079 end
0080 %
0081 %    Temporal spectrum
0082 %
0083 tk = sort(rand(Ni,1))*1000000;
0084 %
0085 ko = ki;emp = zeros(size(ki));
0086 for i = Ni:-1:1,
0087    for k = 1:i-1,
0088       if tk(i) - tk(i-k) < tl+tp,
0089          if tk(i) - tk(i-k) < tl - tp,
0090             a = ki(i-k)/((tp^(alpha+beta))*((beta/alpha)^beta));
0091               t = tp+tk(i)-tk(i-k);
0092               ko(i) = ko(i)+a*(t.^alpha).*((tl-t).^beta);
0093             if tk(i) - tk(i-k) < tp,
0094                a = ki(i)/((tp^(alpha+beta))*((beta/alpha)^beta));
0095                  t = tp-tk(i)+tk(i-k);
0096                  ko(i-k) = ko(i-k)+a*(t.^alpha).*((tl-t).^beta);
0097                if tk(i) - tk(i-k) < tpeak,
0098                   emp(i) = 2;
0099                   emp(i-k) = 1;
0100                else
0101                   emp(i) = 1;
0102                   emp(i-k) = 1;
0103                end
0104             else
0105                emp(i) = 1;
0106             end
0107         else
0108             emp(i) = 1;
0109         end
0110       else
0111          break;
0112       end
0113    end
0114 end
0115 %
0116 [histoko,dummy] = hist_dke_yp(ko,khisto);%    Piled-up spectrum
0117 %
0118 kADC = ko;
0119 if strcmp(tok,'JET'),%JET case only
0120    i = 1;
0121    while i <= Ni,
0122       tt = [];kk = [];
0123       tt(1) = tk(i);
0124       kk(1) = ki(i);
0125       for j = i+1:Ni,
0126          if tk(j) - tk(j-1) < tl,
0127             tt(j-i+1) = tk(j);
0128             kk(j-i+1) = ki(j);
0129          else
0130              break;
0131          end
0132       end
0133       npl = length(kk)-1;
0134       kADC(i:1:i+npl) = jpilupyp(tt,kk,tl,tp,tc,kmin,kmax,nhisto,ko,i);
0135       i = i + npl +1;
0136    end
0137 else %Tore Supra case or other tokamaks
0138    i = 1;
0139    while i <= Ni,
0140       for j = i+1:Ni,
0141          if tk(j)-tk(i)<=tc+tp,
0142                kADC(j) = -1;
0143          else
0144                 break;
0145          end
0146       end
0147       i = j+1;
0148    end
0149 end
0150 %
0151 kkADC = kADC(kADC~=-1);
0152 zz = kADC(find(kADC(:)>0));
0153 [histokADC,dummy] = hist_dke_yp(zz,khisto); %Piled-up + dead time spectrum
0154 %
0155 if pur==0,%Pile-up rejection off
0156    kcode = kkADC;
0157    zz = kcode(find(kcode(:)>0));
0158    [histokcode,dummy] = hist_dke_yp(zz,khisto); % Measured spectrum
0159 elseif pur==1,%Pile-up rejection on
0160    kcode = kADC((kADC~=-1) & (emp~=1));% Measured spectrum
0161    zz = kcode(find(kcode(:)>0));
0162    [histokcode,dummy] = hist_dke_yp(zz,khisto); 
0163 end
0164 
0165 
0166 
0167 
0168 
0169