Context Navigation

V_VSANS_Event_Testing.ipf @ 1046

Last change on this file since 1046 was 1046, checked in by srkline, 6 years ago

changes to properly read and decode VSANS event files, and to process them with a trimmed down (and to be customized) version of the event mode panel.

currently reading and decoding is done without an XOP since the file structure is much simpler than ordela.

TOF mode data can be processed - as this will be needed initially at VSANS to be able to calibrate the wavelength.

File size: 14.9 KB

Line
1	#pragma TextEncoding = "MacRoman"
2	#pragma rtGlobals=3 // Use modern global access method and strict wave access.
3
4
5	//
6	// for the event mode data with the proposed 64 bit structure, I may be able to use Igor for everything.
7	//
8	// Skipping the appropriate bits of the header (after I read them in) may be possible with
9	// either LoadWave (treating the entire wave as 64 bit, unsigned), loading chunks from clipboard?
10	// -- see in LoadWave, the suggestions for "Loading Large Waves"
11	//
12	// or using FBinRead, again, filling a wave, or a chunk of the data, as needed
13	//
14	// Then - since the data has sequential timing, not dependent on rollovers, I can
15	// chunk the data for parallel processing, and piece it back together
16	// after all of the decoding is done.
17	//
18	//
19	// I don't know if it's possible to use a STRUCT definition for each 64 bit word so that I could address
20	// the bytes directly - since that may not work properly in Igor, and I'm not sure how to address the 6 byte section
21	// -possibly with a uchar(6) definition?
22	//
23	//
24	// (from WM help) You can define an array of structures as a field in a structure:
25	// Structure mystruct
26	// STRUCT Point pt[100] // Allowed as a sub-structure
27	// EndStructure
28	//
29	// -- which may be of use to read "blocks" of datq from a file using FBinRead
30	//
31	// -- right now, I'm having issues with being able to "cast" or convert uint64 values to STRUCT
32	// (I don't know how to read from the struct if I can't fill it??)
33	//
34	// TODO:
35	//
36	// (5/2017)
37	// The basic bits of reading work, but will need to be customized to be able to accomodate file names in/out
38	// and especially the number of disabled tubes (although as long as I have the offset, it shouldn't be that
39	// big of an issue.
40	//
41	// -- don't see the struct idea working out. only in real c-code if needed
42	//
43	// -- need to add detector binning to the decoding, to place the counts within the correct panels
44	// -- not sure how this will work with JointHistogram Operation
45	// (split to separate "streams" of values for each detector panel?
46	//
47	//
48	// -- can I efficiently use "sort" (on the tube index) to block the data into groups that can be split
49	// into 4 sets of waves
50	// that can then be binned per panel, using the usual Joint histogram procedures? Works only if the
51	// tube indexing is orderly. if it's a mess, Ill need to try something else (indexed sort?) (replace?)
52	// (manually? ugh.)
53	//
54
55
56	//
57	Structure eventWord
58	uchar eventTime[6]
59	uchar location
60	uchar tube
61	endStructure
62
63
64
65	//
66	//
67	Function V_testBitShift()
68
69	// // /l=64 bit, /U=unsigned
70	// Make/L/U/N=100 eventWave
71	// eventWave = 0
72
73	// for each 64-bit value:
74	// byte 1: tube index [0,191]
75	// byte 2: pixel value [0,127]
76	// bytes 3-8 (= 6 bytes): time stamp in resolution unit
77
78	int64 i64_num,b1,b2,b3,b4,b5,b6,b7,b8
79	int64 i64_ticks,i64_start
80
81	b1=255
82	b3=255
83	b5=255
84	b7=255
85	b2=0
86	b4=0
87	b6=0
88	b8=0
89
90	b7 = b7 << 8
91	b6 = b6 << 16
92	b5 = b5 << 24
93	b4 = b4 << 32
94	b3 = b3 << 40
95	b2 = b2 << 48
96	b1 = b1 << 56
97
98	i64_num = b1+b2+b3+b4+b5+b6+b7+b8
99	printf "%64b\r",i64_num
100
101	return(0)
102	End
103
104	Function V_MakeFakeEvents()
105
106	// // /l=64 bit, /U=unsigned
107	Make/O/L/U/N=10 smallEventWave
108	smallEventWave = 0
109
110	// for each 64-bit value:
111	// byte 1: tube index [0,191]
112	// byte 2: pixel value [0,127]
113	// bytes 3-8 (= 6 bytes): time stamp in resolution unit
114
115	uint64 i64_num,b1,b2,b3,b4,b5,b6,b7,b8
116	uint64 i64_ticks,i64_start
117
118	// b1 = 47
119	// b2 = 123
120	// i64_ticks = 123456789
121	b1 = 41
122	b2 = 66
123	i64_ticks = 15
124
125
126	// b2 = b2 << 48
127	// b1 = b1 << 56
128	//
129	// i64_num = b1+b2+i64_ticks
130
131	// don't shift b1
132	b2 = b2 << 8
133	i64_ticks = i64_ticks << 16
134
135	i64_num = b1+b2+i64_ticks
136
137	printf "%64b\r",i64_num
138	print i64_num
139
140	smallEventWave[0] = i64_num
141
142	return(0)
143	End
144
145	Function V_decodeFakeEvent()
146
147	WAVE w = smallEventWave
148	uint64 val,b1,b2,btime
149	val = w[0]
150
151	// printf "%64b\r",w[0] //wrong (drops the last Å 9 bits)
152	printf "%64b\r",val //correct, assign value to 64bit variable
153	// print w[0] //wrong
154	print val // correct
155
156	// b1 = (val >> 56 ) & 0xFF // = 255, last byte, after shifting
157	// b2 = (val >> 48 ) & 0xFF
158	// btime = val & 0xFFFFFFFFFFFF // = really big number, last 6 bytes
159
160
161	b1 = val & 0xFF
162	b2 = (val >> 8) & 0xFF
163	btime = (val >> 16)
164
165
166	print b1
167	print b2
168	print btime
169
170
171	// //test as struct
172	// Print "as STRUCT"
173	//
174	// STRUCT eventWord s
175	//
176	// s = w[0]
177	//
178	// print s.tube
179	// print s.location
180	// print s.eventTime
181
182
183
184	return(0)
185	End
186
187	//
188	// tested up to num=1e8 successfully
189	//
190	Function V_MakeFakeEventWave(num)
191	Variable num
192
193	Variable ii
194
195
196	// num = 1e3
197
198	// // /l=64 bit, /U=unsigned
199	Make/O/L/U/N=(num) eventWave
200	eventWave = 0
201
202	// for each 64-bit value:
203	// byte 1: tube index [0,191]
204	// byte 2: pixel value [0,127]
205	// bytes 3-8 (= 6 bytes): time stamp in resolution unit
206
207	uint64 i64_num,b1,b2
208	uint64 i64_ticks,i64_start
209
210	i64_start = ticks
211	for(ii=0;ii<num;ii+=1)
212	// sleep/T/C=-1 1 // 6 ticks, approx 0.1 s (without the delay, the loop is too fast)
213	b1 = trunc(abs(enoise(192))) //since truncated, need 192 as highest random to give 191 after trunc
214	b2 = trunc(abs(enoise(128))) // same here, to get results [0,127]
215
216	// i64_ticks = ticks-i64_start
217	i64_ticks = ii+1
218
219	// b2 = b2 << 48
220	// b1 = b1 << 56
221
222	// don't shift b1
223	b2 = b2 << 8
224	i64_ticks = i64_ticks << 16
225
226	i64_num = b1+b2+i64_ticks
227
228	eventWave[ii] = i64_num
229	endfor
230
231
232	return(0)
233	End
234
235
236	//
237	// TODO:
238	// -- can this be multithreaded (eliminating the loop)?
239	//
240	// MultiThread tube = (w >> 56 ) & 0xFF
241	// MultiThread location = (w >> 48 ) & 0xFF
242	// MultiThread eventTime = val & 0xFFFFFFFFFFFF
243	//
244	//
245	//
246	Function V_decodeFakeEventWave(w)
247	Wave w
248
249	s_tic()
250	// WAVE w = eventWave
251	uint64 val,b1,b2,btime
252	val = w[0]
253
254	// printf "%64b\r",w[0] //wrong (drops the last Å 9 bits)
255	// printf "%64b\r",val //correct, assign value to 64bit variable
256	// print w[0] //wrong
257	// print val // correct
258
259	Variable num,ii
260	num=numpnts(w)
261
262	Make/O/L/U/N=(num) eventTime
263	Make/O/U/B/N=(num) tube,location //8 bit unsigned
264
265	for(ii=0;ii<num;ii+=1)
266	val = w[ii]
267
268	// b1 = (val >> 56 ) & 0xFF // = 255, last two bytes, after shifting
269	// b2 = (val >> 48 ) & 0xFF
270	// btime = val & 0xFFFFFFFFFFFF // = really big number, last 6 bytes
271
272	b1 = val & 0xFF
273	b2 = (val >> 8) & 0xFF
274	btime = (val >> 16)
275
276	tube[ii] = b1
277	location[ii] = b2
278	eventTime[ii] = btime
279
280	endfor
281
282	s_toc()
283
284	return(0)
285	End
286
287
288	Function V_writeFakeEventFile(fname)
289	String fname
290
291	WAVE w = eventWave
292	Variable refnum
293
294	String vsansStr="VSANS"
295	Variable revision = 11
296	Variable offset = 26 // no disabled tubes
297	Variable time1 = 2017
298	Variable time2 = 0525
299	Variable time3 = 1122
300	Variable time4 = 3344 // these 4 time pieces are supposed to be 8 bytes total
301	Variable time5 = 3344 // these 5 time pieces are supposed to be 10 bytes total
302	String detStr = "M"
303	Variable volt = 1500
304	Variable resol = 1e7
305
306
307	Open refnum as fname
308
309	FBinWrite refnum, vsansStr
310	FBinWrite/F=2/U refnum, revision
311	FBinWrite/F=2/U refnum, offset
312	FBinWrite/F=2/U refnum, time1
313	FBinWrite/F=2/U refnum, time2
314	FBinWrite/F=2/U refnum, time3
315	FBinWrite/F=2/U refnum, time4
316	FBinWrite/F=2/U refnum, time5
317	FBinWrite refnum, detStr
318	FBinWrite/F=2/U refnum, volt
319	FBinWrite/F=3/U refnum, resol
320
321	FGetPos refnum
322	Print "End of header = ",V_filePos
323	offset = V_filePos
324
325	FSetPos refnum,7
326	FBinWrite/F=2/U refnum, offset //write the correct offset
327
328
329	FSetPos refNum, offset
330
331	FBinWrite refnum, w
332
333	close refnum
334
335	return(0)
336	End
337
338	//
339	// use GBLoadWave to do the reading, then I can do the decoding
340	//
341	Function V_readFakeEventFile(fileName)
342	String filename
343
344	// this reads in uint64 data, to a unit64 wave, skipping 22 bytes
345	// GBLoadWave/B/T={192,192}/W=1/S=22
346	Variable num,refnum
347
348
349	// to read a VSANS event file:
350	//
351	// - get the file name
352	// - read the header (all of it, since I need parts of it) (maybe read as a struct? but I don't know the size!)
353	// - move to EOF and close
354	//
355	// - Use GBLoadWave to read the 64-bit events in
356
357	String vsansStr=""
358	Variable revision
359	Variable offset // no disabled tubes
360	Variable time1
361	Variable time2
362	Variable time3
363	Variable time4 // these 4 time pieces are supposed to be 8 bytes total
364	Variable time5 // these 5 time pieces are supposed to be 10 bytes total
365	String detStr=""
366	Variable volt
367	Variable resol
368
369	vsansStr = PadString(vsansStr,5,0x20) //pad to 5 bytes
370	detStr = PadString(detStr,1,0x20) //pad to 1 byte
371
372	Open/R refnum as filename
373	filename = S_fileName
374
375	s_tic()
376
377	FBinRead refnum, vsansStr
378	FBinRead/F=2/U refnum, revision
379	FBinRead/F=2/U refnum, offset
380	FBinRead/F=2/U refnum, time1
381	FBinRead/F=2/U refnum, time2
382	FBinRead/F=2/U refnum, time3
383	FBinRead/F=2/U refnum, time4
384	FBinRead/F=2/U refnum, time5
385	FBinRead refnum, detStr //NOTE - the example data file Phil sent skipped the detStr (no placeholder!)
386	FBinRead/F=2/U refnum, volt
387	FBinRead/F=3/U refnum, resol
388
389	FStatus refnum
390	FSetPos refnum, V_logEOF
391
392	Close refnum
393
394	// number of data bytes
395	num = V_logEOF-offset
396	Print "Number of data values = ",num/8
397
398	GBLoadWave/B/T={192,192}/W=1/S=(offset) filename // intel, little-endian
399	// GBLoadWave/T={192,192}/W=1/S=(offset) filename // motorola, big-endian
400
401	Duplicate/O $(StringFromList(0,S_waveNames)) V_Events
402	KillWaves/Z $(StringFromList(0,S_waveNames))
403	s_toc()
404
405	Print vsansStr
406	Print revision
407	Print offset
408	Print time1
409	Print time2
410	Print time3
411	Print time4
412	Print time5
413	Print detStr
414	print volt
415	print resol
416
417	return(0)
418	End
419
420	//
421	//
422	//
423	Function V_MakeFakeEventWave_TOF(delayTime,std)
424	Variable delayTime,std
425
426	Variable num,ii,jj,numRepeat
427
428
429	num = 1000
430	numRepeat = 1000
431
432	// delayTime = 50 //microseconds
433	// std = 4 //std deviation, microseconds
434
435	// // /l=64 bit, /U=unsigned
436	Make/O/L/U/N=(num*numRepeat) eventWave
437	eventWave = 0
438
439	Make/O/D/N=(num) arrival
440
441	// for each 64-bit value:
442	// byte 1: tube index [0,191]
443	// byte 2: pixel value [0,127]
444	// bytes 3-8 (= 6 bytes): time stamp in resolution unit
445
446	uint64 i64_num,b1,b2,b3,b4,b5,b6,b7,b8
447	uint64 i64_ticks,i64_start
448
449	// i64_start = ticks
450	i64_ticks = 0
451	for(jj=0;jj<numRepeat;jj+=1)
452	arrival = delayTime + gnoise(std)
453	sort arrival,arrival
454	arrival *= 1000 //milliseconds now
455
456	for(ii=0;ii<num;ii+=1)
457	// sleep/T/C=-1 1 // 6 ticks, approx 0.1 s (without the delay, the loop is too fast)
458	b1 = trunc(abs(enoise(192))) //since truncated, need 192 as highest random to give 191 after trunc
459	b2 = trunc(abs(enoise(128))) // same here, to get results [0,127]
460
461	i64_ticks = trunc(arrival[ii])
462
463	// b2 = b2 << 48
464	// b1 = b1 << 56
465
466	// don't shift b1
467	b2 = b2 << 8
468	i64_ticks = i64_ticks << 16
469
470	i64_num = b1+b2+i64_ticks
471	eventWave[jj*num+ii] = i64_num
472	endfor
473
474	endfor
475
476	return(0)
477	End
478
479
480	// TODO:
481	//
482	// There may be memory issues with this
483	//
484	// -- do I want to do the time binning first?
485	// -- does it really matter?
486	//
487	Function V_SortAndSplitFakeEvents()
488
489	Wave eventTime = root:EventTime
490	Wave location = root:location
491	Wave tube = root:tube
492
493	Sort tube,tube,eventTime,location
494
495	Variable b1,e1,b2,e2,b3,e3,b4,e4
496	FindValue/S=0/I=48 tube
497	b1 = 0
498	e1 = V_Value - 1
499	b2 = V_Value
500	FindValue/S=(b2)/I=96 tube
501	e2 = V_Value - 1
502	b3 = V_Value
503	FindValue/S=(b3)/I=144 tube
504	e3 = V_Value - 1
505	b4 = V_Value
506	e4 = numpnts(tube)-1
507
508	Print b1,e1
509	Print b2,e2
510	Print b3,e3
511	Print b4,e4
512
513	// tube and location become x and y, and can be byte data
514	// eventTime still needs to be 64 bit - when do I convert it to FP?
515	Make/O/B/U/N=(e1-b1+1) tube1,location1
516	Make/O/L/U/N=(e1-b1+1) eventTime1
517
518	Make/O/B/U/N=(e2-b2+1) tube2,location2
519	Make/O/L/U/N=(e2-b2+1) eventTime2
520
521	Make/O/B/U/N=(e3-b3+1) tube3,location3
522	Make/O/L/U/N=(e3-b3+1) eventTime3
523
524	Make/O/B/U/N=(e4-b4+1) tube4,location4
525	Make/O/L/U/N=(e4-b4+1) eventTime4
526
527
528	tube1 = tube[p+b1]
529	tube2 = tube[p+b2]
530	tube3 = tube[p+b3]
531	tube4 = tube[p+b4]
532
533	location1 = location[p+b1]
534	location2 = location[p+b2]
535	location3 = location[p+b3]
536	location4 = location[p+b4]
537
538	eventTime1 = eventTime[p+b1]
539	eventTime2 = eventTime[p+b2]
540	eventTime3 = eventTime[p+b3]
541	eventTime4 = eventTime[p+b4]
542
543
544	KillWaves/Z eventTime,location,tube
545
546	return(0)
547	End
548
549
550
551	// TODO:
552	//
553	// There may be memory issues with this
554	//
555	// -- do I want to do the time binning first?
556	// -- does it really matter?
557	//
558	Function V_SortAndSplitEvents()
559
560
561	SetDataFolder root:Packages:NIST:VSANS:Event:
562
563	Wave eventTime = EventTime
564	Wave location = location
565	Wave tube = tube
566
567	Sort tube,tube,eventTime,location
568
569	Variable b1,e1,b2,e2,b3,e3,b4,e4
570	FindValue/S=0/I=48 tube
571	b1 = 0
572	e1 = V_Value - 1
573	b2 = V_Value
574	FindValue/S=(b2)/I=96 tube
575	e2 = V_Value - 1
576	b3 = V_Value
577	FindValue/S=(b3)/I=144 tube
578	e3 = V_Value - 1
579	b4 = V_Value
580	e4 = numpnts(tube)-1
581
582	Print b1,e1
583	Print b2,e2
584	Print b3,e3
585	Print b4,e4
586
587	// tube and location become x and y, and can be byte data
588	// eventTime still needs to be 64 bit - when do I convert it to FP?
589	Make/O/B/U/N=(e1-b1+1) tube1,location1
590	Make/O/L/U/N=(e1-b1+1) eventTime1
591
592	Make/O/B/U/N=(e2-b2+1) tube2,location2
593	Make/O/L/U/N=(e2-b2+1) eventTime2
594
595	Make/O/B/U/N=(e3-b3+1) tube3,location3
596	Make/O/L/U/N=(e3-b3+1) eventTime3
597
598	Make/O/B/U/N=(e4-b4+1) tube4,location4
599	Make/O/L/U/N=(e4-b4+1) eventTime4
600
601
602	tube1 = tube[p+b1]
603	tube2 = tube[p+b2]
604	tube3 = tube[p+b3]
605	tube4 = tube[p+b4]
606
607	location1 = location[p+b1]
608	location2 = location[p+b2]
609	location3 = location[p+b3]
610	location4 = location[p+b4]
611
612	eventTime1 = eventTime[p+b1]
613	eventTime2 = eventTime[p+b2]
614	eventTime3 = eventTime[p+b3]
615	eventTime4 = eventTime[p+b4]
616
617
618	KillWaves/Z eventTime,location,tube
619
620	return(0)
621	End
622
623
624	//
625	// switch the "active" panel to the selected group (1-4) (5 concatenates them all together)
626	//
627	// copy the set of tubes over to the "active" set that is to be histogrammed
628	// and redimension them to be sure that they are double precision
629	//
630	Function V_SwitchTubeGroup(tubeGroup)
631	Variable tubeGroup
632
633	SetDataFolder root:Packages:NIST:VSANS:Event:
634
635	if(tubeGroup <= 4)
636	Wave tube = $("tube"+num2Str(tubeGroup))
637	Wave location = $("location"+num2Str(tubeGroup))
638	Wave eventTime = $("eventTime"+num2Str(tubeGroup))
639
640	Wave xloc,yLoc,timePt
641
642	KillWaves/Z timePt,xLoc,yLoc
643	Duplicate/O tube xLoc
644	Duplicate/O location yLoc
645	Duplicate/O eventTime timePt
646
647	Redimension/D xLoc,yLoc,timePt
648
649	endif
650
651	if(tubeGroup == 5)
652	Wave xloc,yLoc,timePt
653
654	KillWaves/Z timePt,xLoc,yLoc
655
656	String str = ""
657	str = "tube1;tube2;tube3;tube4;"
658	Concatenate/O/NP str,xloc
659	str = "location1;location2;location3;location4;"
660	Concatenate/O/NP str,yloc
661	str = "eventTime1;eventTime2;eventTime3;eventTime4;"
662	Concatenate/O/NP str,timePt
663
664	Redimension/D xLoc,yLoc,timePt
665	endif
666
667
668	return(0)
669	End
670
671	Proc V_SwitchGroupAndCleanup(num)
672	Variable num
673
674	V_SwitchTubeGroup(num)
675	SetDataFolder root:Packages:NIST:VSANS:Event:
676	Duplicate/O timePt rescaledTime
677	KillWaves/Z OscSortIndex
678	print WaveMax(rescaledTime)
679	root:Packages:NIST:VSANS:Event:gEvent_t_longest = waveMax(rescaledTime)
680
681	SetDataFolder root:
682
683	end
684
685	Function V_count(num)
686	Variable num
687
688	SetDataFolder root:Packages:NIST:VSANS:Event:
689
690	Wave xloc = xloc
691	wave yloc = yloc
692	Variable ii,npt,total=0
693	npt = numpnts(xloc)
694	for(ii=0;ii<npt;ii+=1)
695	if(xloc[ii] == num)
696	total += 1
697	endif
698	if(yloc[ii] == num)
699	total += 1
700	endif
701	endfor
702
703	Print total
704
705	SetDataFolder root:
706
707	end

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format