Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TPCRateLayered.m
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TPCRateLayered.m
1 % bottom up calcualtion of TPC rate using MC approach
2 % Author - Jin Huang <jhuang@bnl.gov>
3 
4 close all
5 
6 clear all
7 
8 myColorMap = colormap;
9 myColorMap(1,:) = [1 1 1];
10 colormap(myColorMap);
11 %% Inputs
12 
13 % n = 1000;
14 % full_rate = 100e3;
15 % trig_rate = 15e3;
16 % trigger_window = 18e-6;
17 
18 % n = 10000000;
19 n = 50000;
20 % n = 500000;
21 % full_rate = 250e3;
22 full_rate = 100e3;
23 % full_rate = 150e3;
24 % full_rate = 100e3;
25 % full_rate = 50e3;
26 trig_rate = 15e3;
27 % trigger_window = 4e-6;
28 trigger_window = 13e-6;
29 % trigger_window = 17.5e-6;
30 % trigger_window = 35e-6;
31 
32 TargetEta = 1.1;
33 % TargetEta = 3;
34 
35 nRing = 40;
36 % nRing = 48;
37 
38 SaveName = sprintf('TPCRateLayered_%.0fHzCol_%.0fHzTrig_%.0fusDrift_nRing%d',full_rate,trig_rate,trigger_window*1e6, nRing);
39 
40 %% Generic constants
41 
42 
43 minR = 30;
44 maxR = 78;
45 maxZ = 105;
46 dNdeta = 2 * 180 * 2; % Pre-CDR table 3.3 with effective factor x2 and two signs x2
47 bitPerHit =3*5*10 * 1.4; % Pre-CDR table 3.3 with effective factor
48 DAMCompressionFactor = 1.02 * 0.5 * 0.6; % Repacking, clustering, compression
49 
50 BCO = 10e6;
51 DataPlotRangeBCO =10000;
52 
53 TimeSpan= n/full_rate;
54 TimeSpanBCO = int64(TimeSpan*BCO);
55 
56 
57 TriggerBCO = poissrnd(trig_rate/BCO, TimeSpanBCO, 1)';
58 CollisionBCO =poissrnd((full_rate - trig_rate)/BCO ,TimeSpanBCO, 1)' + TriggerBCO;
59 
60 
61 %% Per triggger
62 
63 PerTriggerBCO = int64(trigger_window * BCO);
64 
65 RRing = repmat(linspace(minR, maxR, nRing )',1,PerTriggerBCO);
66 zBCO = repmat(linspace(maxZ, 0, PerTriggerBCO),nRing,1);
67 
68 MinEtaRingBCO = atanh( (zBCO - 10) ./ sqrt((zBCO - 10).^2 + RRing.^2) );
69 % MinEtaRingBCO = atanh( (zBCO - 0) ./ sqrt((zBCO - 0).^2 + RRing.^2) );
70 dEtaRingBCO = atanh( (zBCO + maxZ/double(PerTriggerBCO) ) ./ sqrt((zBCO + maxZ/double(PerTriggerBCO) ).^2 + RRing.^2) ) - atanh( zBCO ./ sqrt(zBCO.^2 + RRing.^2) );
71 dataBitRingBCO = bitPerHit .* dNdeta .* dEtaRingBCO; % Pre-CDR table 3.3 with effective factor
72 
73 TriggerMask = double(MinEtaRingBCO < TargetEta);
74 dataBitRingBCOAfterMask = dataBitRingBCO.*TriggerMask;
75 
76 figure('name','PerEventData','PaperPositionMode','auto', ...
77  'position',[100,0,2400,800]) ;
78 
79 subplot(1,2,1)
80 
81 
82 imagesc(zBCO(1,:), RRing(:,1), dataBitRingBCO);
83 c = colorbar;
84 
85 hold on;
86 plot([10 maxZ+10], [0 tan(acos(tanh(TargetEta))) * maxZ],'k--');
87 
88 set(gca,'YDir','normal');
89 set(gca,'YLim',[0,80]);
90 set(gca,'XLim',[0,110])
91 box on
92 grid on
93 daspect([1 1 1])
94 
95 xlabel('|z| position (cm)','FontSize',14);
96 ylabel('R position (cm)','FontSize',14);
97 c.Label.String = 'FEE data (bit) per MB collision per |z|-R bin (integrated over z-sign and azimuth)';
98 c.Label.FontSize = 14;
99 
100 title(sprintf('FEE data (bit) per MB collision per |z|-R bin, total = %.1f Mbit', sum(sum(dataBitRingBCO))/1e6),'FontSize',16);
101 
102 
103 subplot(1,2,2)
104 
105 
106 imagesc(zBCO(1,:), RRing(:,1), dataBitRingBCOAfterMask);
107 c = colorbar;
108 hold on;
109 plot([10 maxZ+10], [0 tan(acos(tanh(TargetEta))) * maxZ],'k--');
110 
111 set(gca,'YDir','normal');
112 set(gca,'YLim',[0,80]);
113 set(gca,'XLim',[0,110])
114 box on
115 grid on
116 daspect([1 1 1])
117 
118 xlabel('|z| position (cm)','FontSize',14);
119 ylabel('R position','FontSize',14);
120 c.Label.String = 'FEE data (bit) per MB collision per |z|-R bin (integrated over z-sign and azimuth)';
121 c.Label.FontSize = 14;
122 
123 title(sprintf('FEE data (bit) after DAM acceptance filtering, total = %.1f Mbit', sum(sum(dataBitRingBCOAfterMask))/1e6),'FontSize',16);
124 
125 
126 
127 
128 SaveCavas('TPCRateLayered',gcf);
129 
130 
131 %% Run the DAQ
132 
133 DataBCO = conv2(CollisionBCO,dataBitRingBCO);
134 TriggerAcceptBCO = conv2(TriggerBCO,TriggerMask);
135 ThrottleAcceptBCO = double(TriggerAcceptBCO>0);
136 
137 PlotDataBCO = DataBCO(:,1:DataPlotRangeBCO);
138 PlotThrottleDataBCO = DataBCO(:,1:DataPlotRangeBCO).*ThrottleAcceptBCO(:,1:DataPlotRangeBCO);
139 
140 TotalDataRate = sum(sum(DataBCO))/TimeSpan;
141 TriggerDataRate = sum(sum(DataBCO.*TriggerAcceptBCO))/TimeSpan;
142 ThrottleDataRate = sum(sum(DataBCO.*ThrottleAcceptBCO))/TimeSpan;
143 
144 %% Event stream display
145 
146 figure('name','DataStream','PaperPositionMode','auto', ...
147  'position',[100,0,2400,800]) ;
148 
149 subplot(2,1,1)
150 
151 imagesc( DataBCO(:,1:DataPlotRangeBCO) );
152 c = colorbar;
153 hold on;
154 set(gca,'YDir','normal');
155 set(gca,'YLim',[-4,nRing+1]);
156 
157 xlabel('BCO bin (100 ns)','FontSize',14);
158 ylabel('Radial layer number','FontSize',14);
159 c.Label.String = 'FEE data (bit)';
160 c.Label.FontSize = 14;
161 title(sprintf('FEE data input to DAM. Rate = %.0f Gbps @ %.0f kHz Collision, %.0f kHz Trigger %.0f us Drift',...
162  sum(sum(DataBCO))/TimeSpan/1e9,full_rate/1e3,trig_rate/1e3,trigger_window*1e6),'FontSize',16);
163 
164 for i = 1:DataPlotRangeBCO
165 
166  if (TriggerBCO(i) >0)
167 
168  plot([i i], [-3 nRing],'r--','LineWidth',2);
169  plot([i,i+PerTriggerBCO], [-3 -3], 'r-', 'LineWidth', 2)
170  plot([i,i+PerTriggerBCO], [-1 -1], 'k-', 'LineWidth', 2)
171 
172  elseif (CollisionBCO(i) >0)
173 
174  plot([i i], [-1 nRing],'k--','LineWidth',2);
175  plot([i,i+PerTriggerBCO], [-1 -1], 'k-', 'LineWidth', 2)
176  end
177 
178 
179 end
180 
181 
182 
183 subplot(2,1,2)
184 
185 imagesc( PlotThrottleDataBCO );
186 c = colorbar;
187 hold on;
188 set(gca,'YDir','normal');
189 set(gca,'YLim',[-4,nRing+1]);
190 
191 xlabel('BCO bin (100 ns)','FontSize',14);
192 ylabel('Radial layer number','FontSize',14);
193 c.Label.String = 'FEE data equavelent (bit)';
194 c.Label.FontSize = 14;
195 title(sprintf('DAM data output: Throttled data rate = %.0f Gbps, Triggered data rate = %.0f Gbps ',...
196  ThrottleDataRate * DAMCompressionFactor/1e9 ,...
197  TriggerDataRate * DAMCompressionFactor/1e9...
198  ),'FontSize',16);
199 
200 for i = 1:DataPlotRangeBCO
201 
202  if (TriggerBCO(i) >0)
203 
204  plot([i i], [-3 nRing],'r--','LineWidth',2);
205  plot([i,i+PerTriggerBCO], [-3 -3], 'r-', 'LineWidth', 2)
206  plot([i,i+PerTriggerBCO], [-1 -1], 'k-', 'LineWidth', 2)
207 
208  elseif (CollisionBCO(i) >0)
209 
210  plot([i i], [-1 nRing],'k--','LineWidth',2);
211  plot([i,i+PerTriggerBCO], [-1 -1], 'k-', 'LineWidth', 2)
212  end
213 
214 
215 end
216 
217 SaveCavas(SaveName,gcf);
218 
219 
220 %% Event stream display - 1D
221 
222 figure('name','DataStream1D','PaperPositionMode','auto', ...
223  'position',[100,0,2400,800]) ;
224 
225 subplot(2,1,1)
226 
227 bar( PlotDataBCO(nRing,:) ,1);
228 % c = colorbar;
229 hold on;
230 set(gca,'YDir','normal');
231 set(gca, 'XLim',[0 DataPlotRangeBCO]);
232 ylim = get(gca,'YLim');
233 y_max = ylim(2);
234 
235 xlabel('BCO bin (100 ns)','FontSize',14);
236 ylabel('Data in last TPC layer (bit/BCO)','FontSize',14);
237 % c.Label.String = 'FEE data (bit)';
238 % c.Label.FontSize = 14;
239 title(sprintf('FEE data input to DAM. Rate = %.0f Gbps @ %.0f kHz Collision, %.0f kHz Trigger %.0f us Drift',...
240  sum(sum(DataBCO))/TimeSpan/1e9,full_rate/1e3,trig_rate/1e3,trigger_window*1e6),'FontSize',16);
241 
242 for i = 1:DataPlotRangeBCO
243 
244  if (TriggerBCO(i) >0)
245 
246  plot([i i], [-.15 1].*y_max,'r--','LineWidth',1);
247  plot([i,i+PerTriggerBCO], [-.15 -.15].*y_max, 'r-', 'LineWidth', 1)
248  plot([i,i+PerTriggerBCO], [-.1 -.1].*y_max, 'k-', 'LineWidth', 1)
249 
250  elseif (CollisionBCO(i) >0)
251 
252  plot([i i], [-.1 1].*y_max,'k--','LineWidth',1);
253  plot([i,i+PerTriggerBCO], [-.1 -.1].*y_max, 'k-', 'LineWidth', 1)
254  end
255 
256 
257 end
258 
259 
260 
261 subplot(2,1,2)
262 
263 bar( PlotThrottleDataBCO(nRing,:), 1 );
264 % c = colorbar;
265 hold on;
266 set(gca,'YDir','normal');
267 set(gca, 'XLim',[0 DataPlotRangeBCO]);
268 % set(gca,'YLim',[-4,nRing+1]);
269 
270 xlabel('BCO bin (100 ns)','FontSize',14);
271 ylabel('Data in last TPC layer (bit/BCO)','FontSize',14);
272 % c.Label.String = 'FEE data equavelent (bit)';
273 % c.Label.FontSize = 14;
274 title(sprintf('DAM data output: Throttled data rate = %.0f Gbps, Triggered data rate = %.0f Gbps ',...
275  ThrottleDataRate * DAMCompressionFactor/1e9 ,...
276  TriggerDataRate * DAMCompressionFactor/1e9...
277  ),'FontSize',16);
278 
279 for i = 1:DataPlotRangeBCO
280 
281  if (TriggerBCO(i) >0)
282 
283  plot([i i], [-.15 1].*y_max,'r--','LineWidth',1);
284  plot([i,i+PerTriggerBCO], [-.15 -.15].*y_max, 'r-', 'LineWidth', 1)
285  plot([i,i+PerTriggerBCO], [-.1 -.1].*y_max, 'k-', 'LineWidth', 1)
286 
287  elseif (CollisionBCO(i) >0)
288 
289  plot([i i], [-.1 1].*y_max,'k--','LineWidth',1);
290  plot([i,i+PerTriggerBCO], [-.1 -.1].*y_max, 'k-', 'LineWidth', 1)
291  end
292 
293 
294 end
295 
296 
297 SaveCavas(SaveName,gcf);
298 
299 %% FEE Data Rate
300 
301 DataRadialLayerRate = sum(DataBCO, 2) / TimeSpan;
302 
303 if nRing == 40
304 
305  FEEDataRate = [sum( DataRadialLayerRate(1:8))/5;
306  sum( DataRadialLayerRate((8+1):(8+16)))/8;
307  sum( DataRadialLayerRate((8+1+16):(8+16+16)))/12]...
308  /12/2;%12 sectors and two sides
309 elseif nRing == 48
310 
311  FEEDataRate = [sum( DataRadialLayerRate(1:16))/6;
312  sum( DataRadialLayerRate((16+1):(16+16)))/8;
313  sum( DataRadialLayerRate((16+1+16):(16+16+16)))/12]...
314  /12/2;%12 sectors and two sides
315 
316 else
317  assert(0);
318 end
319 
320 
321 %% FEE Data Rate Plot
322 
323 figure('name','FEEDataRate','PaperPositionMode','auto', ...
324  'position',[100,0,800,400]) ;
325 
326 
327 bar( FEEDataRate/1e9 );
328 % c = colorbar;
329 hold on;
330 set(gca,'YDir','normal');
331 % set(gca, 'XLim',[0 DataPlotRangeBCO]);
332 ylim = get(gca,'YLim');
333 y_max = ylim(2);
334 
335 xlabel('Module Number','FontSize',16);
336 ylabel('Average data rate per FEE (Gbps)','FontSize',16);
337 % c.Label.String = 'FEE data (bit)';
338 % c.Label.FontSize = 14;
339 title(sprintf('FEE data input to DAM. Total = %.0f Gbps @ %.0f kHz Collision, %.0f kHz Trigger %.0f us Drift',...
340  sum(sum(DataBCO))/TimeSpan/1e9,full_rate/1e3,trig_rate/1e3,trigger_window*1e6),'FontSize',12);
341 
342 SaveCavas(SaveName,gcf);
343 
344 %% Collision pile up histogram
345 
346 nPileUpTrig = zeros( sum(TriggerBCO(PerTriggerBCO+1:(TimeSpanBCO-PerTriggerBCO))) ,1) ;
347 nPileUpRnd = zeros( TimeSpanBCO - 2*PerTriggerBCO ,1) ;
348 
349 iTrig = 1;
350 for i = PerTriggerBCO+1:(TimeSpanBCO-PerTriggerBCO)
351 
352  if (TriggerBCO(i));
353 
354  nPileUpTrig(iTrig) = sum( CollisionBCO((i-PerTriggerBCO):(i+PerTriggerBCO)) );
355 
356  iTrig = iTrig +1;
357  end
358 
359  nPileUpRnd(i- PerTriggerBCO + 1) = sum( CollisionBCO((i-PerTriggerBCO):(i+PerTriggerBCO)) );
360 
361 
362 end
363 %%
364 
365 figure('name','nPileUp','PaperPositionMode','auto', ...
366  'position',[100,0,2000,800]) ;
367 
368 subplot(1,2,1)
369 
370 hist(nPileUpTrig,0:20);
371 
372 set(gca,'XLim',[-1,11]);
373 xlabel('Number of collisions in TPC drift window','FontSize',14);
374 ylabel('Count per bin','FontSize',14);
375 title(sprintf('Collision trigger: <# collision in TPC window> = %.2f @ %.0f kHz Collision, %.0f us Drift',...
376  mean(nPileUpTrig),...
377  full_rate/1e3,trigger_window*1e6),'FontSize',16);
378 
379 subplot(1,2,2)
380 
381 hist(nPileUpRnd,0:20);
382 
383 set(gca,'XLim',[-1,11]);
384 xlabel('Number of collisions in TPC drift window','FontSize',14);
385 ylabel('Count per bin','FontSize',14);
386 title(sprintf('Random trigger: <# collision in TPC window> = %.2f @ %.0f kHz Collision, %.0f us Drift',...
387  mean(nPileUpRnd),...
388  full_rate/1e3,trigger_window*1e6),'FontSize',16);
389 
390 SaveCavas(SaveName,gcf);
391 
392 %%
393 
394 % fprintf('Throttled event / total = %.3f; Throttled data / total = %.3f; Triggered event / total = %.3f; Triggered data / total = = %.3f\n',...
395 % n_recorded/n, recorded_data/n, n_triggered./n, triggered_data/n);
396 fprintf('------------\n');
397 
398 fprintf('full_rate = %.0f kHz; trig_rate= %.0f kHz; trigger_window= %.1f us; Time span = %.3f s \n',...
399  full_rate/1e3,trig_rate/1e3, trigger_window*1e6, TimeSpan );
400 
401 fprintf('per event effective track = %.0f , per event FEE data = %.0f bits, total FEE data rate = %.0f Gbps \n',...
402  mean(sum(dEtaRingBCO,2)*180*2), sum(sum(dataBitRingBCO)), sum(sum(DataBCO))/TimeSpan/1e9 );
403 
404 fprintf('Trigger rate*drift window = %.2f;Full rate*drift window= %.2f;Trigger rate/full rate = %.2f (input)/%.2f(MC); \n',...
405  trig_rate*trigger_window,full_rate*trigger_window,trig_rate/full_rate, sum(TriggerBCO)/sum(CollisionBCO) );
406 
407 fprintf('throttled data / total = %.3f; Triggered data / total = %.3f; throttled/trigger = %.3f; Triggered data/Per event data = %.3f \n',...
408  ThrottleDataRate/TotalDataRate,...
409  TriggerDataRate/TotalDataRate,...
410  ThrottleDataRate/TriggerDataRate,...
411  sum(sum(DataBCO.*TriggerAcceptBCO)) ./ (sum(sum(dataBitRingBCOAfterMask)).*sum(TriggerBCO))...
412  );
413 % sum(sum(DataBCO.*TriggerAcceptBCO)) ./ (sum(sum(dataBitRingBCOAfterMask)).*sum(CollisionBCO))...
414 
415 fprintf('throttled data rate = %.0f Gbps; Triggered data rate = %.0f Gbps\n',...
416  ThrottleDataRate * DAMCompressionFactor/1e9 ,...
417  TriggerDataRate * DAMCompressionFactor/1e9);
418 fprintf('------------\n');
419 
420 % fprintf('throttled event / total = %.3f; Triggered event / total = = %.3f; throttled/trigger = %.3f\n',...
421 % recorded_data/n, triggered_data/n, recorded_data/triggered_data );
422