0001 function [Ni_s,khisto1,Tphi_s,eTphi_s,Qi_s,Tpho_s,eTpho_s,Qo_s,TphADC_s,eTphADC_s,QADC_s,Tphcode_s,eTphcode_s,Qcode_s] = pileup_dke_yp(pileupsim,diag)
0002
0003
0004
0005
0006
0007
0008 if nargin < 2,
0009 error('Not enough input arguments for pileup_dke_yp');
0010 end
0011
0012 pur = pileupsim.pur;
0013 nreset = pileupsim.nreset;
0014 typspec = pileupsim.typspec;
0015 Ni0 = pileupsim.Ni0;
0016 ech = pileupsim.ech;
0017 k0 = pileupsim.k0;
0018 modeMC = pileupsim.modeMC;
0019 NiMC = pileupsim.NiMC;
0020 seuil = pileupsim.seuil;
0021 p_opt = pileupsim.p_opt;
0022 filename = pileupsim.filename;
0023
0024 tok = diag.tok;
0025 namediag = diag.namediag;
0026 theta = diag.theta;
0027 tc0 = diag.tc0;
0028 tpeak = diag.tpeak;
0029 tl0 = diag.tl0;
0030 tp0 = diag.tp0;
0031 beta = diag.beta;
0032 kmin = diag.kmin;
0033 kmin0 = diag.kmin0;
0034 kmax = diag.kmax;
0035 nhisto = diag.nhisto;
0036
0037 Ni = [1:500:10*Ni0+1];
0038
0039 tp = tp0*theta;
0040 tl = tl0*theta;
0041 tc = tc0;
0042 if tc < tl-tp, tc1 = tl-tp;else tc1 = tc;end
0043
0044 Ncc = Ni./(exp((tp+tl)*Ni*1e-6)+Ni*(tc1-(tl-tp))*1e-6);
0045 Nc = Ni./exp((tp+tl)*Ni*1e-6);
0046
0047 delta = (Ni-Ncc)./Ni;
0048 index = find(delta<=seuil);
0049
0050 figure('Name','Pile-up level'),
0051 [ax] = graph1D_jd(Ni,Ncc,0,0,'Input count rate (photon/s)','Output count rate (cps/s)',['Pile-up effect ',tok,' ',namediag],NaN,[0,max(Ni)*1.01],[0,max(Nc)*1.5],'-','none','r',2);
0052 [ax] = graph1D_jd(Ni,Nc,0,0,'Input count rate (photon/s)','Output count rate (cps/s)',['Pile-up effect ',tok,' ',namediag],NaN,[0,max(Ni)*1.01],[0,max(Nc)*1.5],'--','none','g',2);
0053 [ax] = graph1D_jd(Ni,Ni,0,0,'Input count rate (photon/s)','Output count rate (cps/s)',['Pile-up effect ',tok,' ',namediag],NaN,[0,max(Ni)*1.01],[0,max(Nc)*1.5],'-.','none','k',2);
0054 legend('Pile-up + ADC','Pile-up only','Ideal','Location','South');
0055 [ax] = graph1D_jd([Ni(max(index)),Ni(max(index))],[1,max(Nc)*1.5],0,0,'Input count rate (pulse/s)','Output count rate (cps/s)',['Pile-up effect ',tok,' ',namediag],NaN,[0,max(Ni)*1.01],[0,max(Nc)*1.5],'--','none','b',1);
0056 [ax] = graph1D_jd([1,max(Ni)],[Ncc(max(index)),Ncc(max(index))],0,0,'Input count rate (pulse/s)','Output count rate (cps/s)',['Pile-up effect ',tok,' ',namediag],NaN,[0,max(Ni)*1.01],[0,max(Nc)*1.5],'--','none','b',1);
0057 ht = text(0,1.05*Ncc(max(index)),[' ',int2str(pileupsim.seuil*100),'% pile-up']);set(ht,'fontsize',20);
0058 axis('square');
0059 print_jd(p_opt,[filename,'_CountRate']);
0060
0061 disp('----------------------------------------------------------------------');
0062 info_dke_yp(1,'Pile-up at high count rate');
0063 disp('----------------------------------------------------------------------');
0064 disp([' - Shaping time (theta) = ',num2str(theta),' 1e-6 (s)']);
0065 disp([' - tp (peak time, approx. 2*theta) = ',num2str(tp),' 1e-6 (s)']);
0066 disp([' - tl (pulse length, approx. 6*theta) = ',num2str(tl),' 1e-6 (s)']);
0067 disp([' - tc (coding dead-time) = ',num2str(tc),' 1e-6 (s)']);
0068 disp([' - Maximum theoretical output count rate without coding dead-time (Nc) = ',num2str(max(Nc)),' (cps/s)']);
0069 disp([' - Maximum theoretical output count rate with coding dead-time(Ncc) = ',num2str(max(Ncc)),' (cps/s)']);
0070 disp([' - For ',num2str(seuil*100),' (%) pile-up: ']);
0071 disp([' - Ni = ',int2str(Ni(max(index))),' (cps/s)']);
0072 disp([' - No (Unpiled-up) = ',int2str(Ncc(max(index))),' (cps/s)']);
0073
0074
0075
0076 disp('----------------------------------------------------------------------');
0077 info_dke_yp(1,['Monte-Carlo study of pile-up and dead time for ',tok,'-',namediag])
0078 disp('----------------------------------------------------------------------');
0079
0080
0081
0082 alpha = beta*tp/(tl-tp);
0083 a = k0/((tp^(alpha+beta))*((beta/alpha)^beta));
0084 t = [0:tl/100:tl];y = a*(t.^alpha).*((tl-t).^beta);
0085 figure('Name','Amplifier pulse shape'),
0086 [ax] = graph1D_jd(t,y,0,0,'Time (microsecond)','E (keV)',['Amplifier pulse shape ',tok,' ',namediag],NaN,'','','-','none','r',2);
0087 print_jd(p_opt,[filename,'_PulseModel']);
0088
0089
0090
0091 Ni_in = Ni(max(index));
0092 [ki,ko,kADC,kkADC,emp,kcode,khisto,khisto1,histoki,histoko,histokADC,histokcode] = fpilup_dke_yp(diag,nreset,Ni_in,kmin0,k0,ech,typspec,pur);
0093 khisto1
0094 khisto
0095
0096
0097
0098 maxfig = 10^(round(log(max(hist_dke_yp(ki,khisto)))/log(10))+1);
0099
0100 figure('Name','Pile-up effect/Spectrum IN'),
0101 hist_dke_yp(ki,khisto,'Pulse energy E (keV)','',['Spectrum IN, N_{tot} = ',int2str(length(ki)),' (cps/s)'],[0,kmax],[1,maxfig],'m','-',2);
0102 if typspec==2,
0103 kint = (khisto1>=diag.khistotphmin & khisto1 <=diag.khistotphmax);
0104 [Tphi,eTphi,Ai,eAi,Qi] = tphfit_dke_yp(khisto1(kint)',histoki(kint)',sqrt(histoki(kint)'));
0105 pfit = polyfit(khisto1,log(histoki),1);
0106 krefi = -1/pfit(1);afiti = exp(pfit(2));
0107 hold on
0108 [ax] = graph1D_jd(khisto1,afiti*exp(-khisto1/krefi),0,1,'k (keV)','','Spectrum IN',NaN,[0,kmax],[1,maxfig],'--','none','m',1);
0109 hold off;
0110 legend('Monte-Carlo','Fit')
0111 end
0112 print_jd(p_opt,[filename,'_IN_PulseSpectrum']);
0113
0114
0115
0116 figure('Name','Pile-up effect/Spectrum OUT'),
0117 hist_dke_yp(ko,khisto,'Pulse energy E (keV)','',['Spectrum OUT, N_{tot} = ',int2str(length(ko)),' (cps/s)'],[0,kmax],[1,maxfig],'m','-',2);
0118 if typspec==2,
0119 [Tpho,eTpho,Ao,eAo,Qo] = tphfit_dke_yp(khisto1(kint)',histoko(kint)',sqrt(histoko(kint)'));
0120 pfit = polyfit(khisto1,log(histoko),1);
0121 krefo = -1/pfit(1);afito = exp(pfit(2));
0122 hold on
0123 [ax] = graph1D_jd(khisto1,afito*exp(-khisto1/krefo),0,1,'k (keV)','','Spectrum OUT',NaN,[0,kmax],[1,maxfig],'--','none','m',1);
0124 hold off;
0125 legend('Monte-Carlo','Fit')
0126 end
0127 print_jd(p_opt,[filename,'_OUT_PulseSpectrum']);
0128
0129
0130
0131 figure('Name','Pile-up effect/Spectrum OUT + coding dead-time'),
0132 hist_dke_yp(kkADC,khisto,'Pulse energy E (keV)','',['Spectrum OUT + coding dead-time, N_{tot} = ',int2str(length(ki)),' (cps/s)'],[0,kmax],[1,maxfig],'m','-',2);
0133 if typspec==2,
0134 [TphADC,eTphADC,AADC,eAADC,QADC] = tphfit_dke_yp(khisto1(kint)',histokADC(kint)',sqrt(histokADC(kint)'));
0135 pfit = polyfit(khisto1,log(histokADC),1);
0136 krefADC = -1/pfit(1);afitADC = exp(pfit(2));
0137 hold on
0138 [ax] = graph1D_jd(khisto1,afitADC*exp(-khisto1/krefADC),0,1,'k (keV)','','Spectrum OUT+ADC',NaN,[0,kmax],[1,maxfig],'--','none','m',1);
0139 hold off;
0140 legend('Monte-Carlo','Fit')
0141 end
0142 print_jd(p_opt,[filename,'_OUTWithDeadTime_PulseSpectrum']);
0143
0144
0145
0146 if pur==0,
0147 titlemode = ['Measured spectrum (PUR OFF), N_{tot} = ',int2str(length(kcode)),' (cps/s)'];
0148 elseif pur==1,
0149 titlemode = ['Measured spectrum (PUR ON), N_{tot} = ',int2str(length(kcode)),' (cps/s)'];
0150 end
0151
0152 figure('Name','Pile-up effect/Measured spectrum'),
0153 hist_dke_yp(kcode,khisto,'Pulse energy E (keV)','',titlemode,[0,kmax],[1,maxfig],'m','-',2);
0154 if typspec==2,
0155 [Tphcode,eTphcode,Acode,eAcode,Qcode] = tphfit_dke_yp(khisto1(kint)',histokcode(kint)',sqrt(histokcode(kint)'));
0156 pfit = polyfit(khisto1,log(histokcode),1);
0157 krefcode = -1/pfit(1);afitcode = exp(pfit(2));
0158 hold on
0159 [ax] = graph1D_jd(khisto1,afitcode*exp(-khisto1/krefcode),0,1,'k (keV)','',titlemode,NaN,[0,kmax],[1,maxfig],'--','none','m',1);
0160 hold off;
0161 legend('Monte-Carlo','Fit')
0162 end
0163 print_jd(p_opt,[filename,'_MeasuredSpectrum']);
0164
0165
0166
0167 info_dke_yp(2,['Number of incoming pulses : ',int2str(Ni_in),' (pulses/s)'])
0168 info_dke_yp(2,['Fraction of piled pulses : ',num2str(100*(sum(emp==1)+sum((emp==2)/2))*ech/Ni_in),' (%)'])
0169 info_dke_yp(2,['Fraction of lost pulses (coding dead-time) : ',num2str(100*sum(kADC==-1)*ech/Ni_in),' (%)'])
0170 skcode = size(kcode);
0171
0172 if pur==0,
0173 info_dke_yp(2,['Fraction of encoded pulses (PUR OFF) : ',num2str(100*skcode(1)*ech/Ni_in),' (%)'])
0174 elseif pur==1,
0175 info_dke_yp(2,['Fraction of encoded pulses (PUR ON) : ',num2str(100*skcode(1)*ech/Ni_in),' (%)'])
0176 end
0177
0178 if typspec==2,
0179 info_dke_yp(2,'Pile-up effect on the measured energy spectrum')
0180 disp([' - Photon energy interval for Tph calculations: ',num2str(min(khisto1(kint))),' - ',num2str(max(khisto1(kint))),' (keV)']);
0181 disp([' - Reference energy k0 (exponential photon spectrum) = ',num2str(k0),' (keV)']);
0182 disp([' - Tph (IN) = ',num2str(Tphi),' +/- ',num2str(eTphi),' (keV), confidence level: ',num2str(Qi)]);
0183 disp([' - Tph (OUT) = ',num2str(Tpho),' +/- ',num2str(eTpho),' (keV), confidence level: ',num2str(Qo)]);
0184 disp([' - Tph (OUT + coding dead-time) = ',num2str(TphADC),' +/- ',num2str(eTphADC),' (keV), confidence level: ',num2str(QADC)]);
0185 if pur==1,
0186 disp([' - Tph (CODED, PUR ON)= ',num2str(Tphcode),' +/- ',num2str(eTphcode),' (keV), confidence level: ',num2str(Qcode)]);
0187 elseif pur==0,
0188 disp([' - Tph (CODED, PUR OFF)= ',num2str(Tphcode),' +/- ',num2str(eTphcode),' (keV), confidence level: ',num2str(Qcode)]);
0189 end
0190 end
0191
0192 Ni_s(1) = Ni_in;
0193 Tphi_s(1) = Tphi;
0194 eTphi_s(1) = eTphi;
0195 Qi_s(1) = Qi;
0196 Tpho_s(1) = Tpho;
0197 eTpho_s(1) = eTpho;
0198 Qo_s(1) = Qo;
0199 TphADC_s(1) = TphADC;
0200 eTphADC_s(1) = eTphADC;
0201 QADC_s(1) = QADC;
0202 Tphcode_s(1) = Tphcode;
0203 eTphcode_s(1) = eTphcode;
0204 Qcode_s(1) = Qcode;
0205
0206 if modeMC == 2,
0207 for ii = 1:length(NiMC),
0208 [ki,ko,kADC,kkADC,emp,kcode,khisto,khisto1,histoki,histoko,histokADC,histokcode] = fpilup_dke_yp(diag,nreset,NiMC(ii),kmin0,k0,ech,typspec,pur);
0209 if typspec==2,
0210 kint = (khisto1>=diag.khistotphmin & khisto1 <=diag.khistotphmax);
0211 [Tphi,eTphi,Ai,eAi,Qi] = tphfit_dke_yp(khisto1(kint)',histoki(kint)',sqrt(histoki(kint)'));
0212 [Tpho,eTpho,Ao,eAo,Qo] = tphfit_dke_yp(khisto1(kint)',histoko(kint)',sqrt(histoko(kint)'));
0213 [TphADC,eTphADC,AADC,eAADC,QADC] = tphfit_dke_yp(khisto1(kint)',histokADC(kint)',sqrt(histokADC(kint)'));
0214 [Tphcode,eTphcode,Acode,eAcode,Qcode] = tphfit_dke_yp(khisto1(kint)',histokcode(kint)',sqrt(histokcode(kint)'));
0215
0216 Ni_s(ii+1) = NiMC(ii);
0217 Tphi_s(ii+1) = Tphi;
0218 eTphi_s(ii+1) = eTphi;
0219 Qi_s(ii+1) = Qi;
0220 Tpho_s(ii+1) = Tpho;
0221 eTpho_s(ii+1) = eTpho;
0222 Qo_s(ii+1) = Qo;
0223 TphADC_s(ii+1) = TphADC;
0224 eTphADC_s(ii+1) = eTphADC;
0225 QADC_s(ii+1) = QADC;
0226 Tphcode_s(ii+1) = Tphcode;
0227 eTphcode_s(ii+1) = eTphcode;
0228 Qcode_s(ii+1) = Qcode;
0229 info_dke_yp(2,['Monte-Carlo calculations for a count rate of ',int2str(NiMC(ii)),' ph/s done.']);
0230 end
0231 end
0232 end
0233