source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/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//
57Structure eventWord
58        uchar eventTime[6]
59        uchar location
60        uchar tube
61endStructure
62
63
64
65//
66//
67Function 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)
102End
103
104Function 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)
143End
144
145Function 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)
185End
186
187//
188// tested up to num=1e8 successfully
189//
190Function 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)
233End
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//
246Function V_decodeFakeEventWave(w)
247        Wave w
248
249s_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
282s_toc()
283               
284        return(0)
285End
286
287
288Function 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)
336End
337
338//
339// use GBLoadWave to do the reading, then I can do the decoding
340//
341Function 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
375s_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))
403s_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)
418End
419
420//
421//
422//
423Function 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)
477End
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//
487Function 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)
547End
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//
558Function 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)
621End
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//
630Function 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)
669End
670
671Proc 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
683end
684
685Function 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
707end
Note: See TracBrowser for help on using the repository browser.