AlcapDAQ  1
show_pulses.cxx
Go to the documentation of this file.
1 #define _FILE_OFFSET_BITS 64
2 
3 /*
4  * slow_extract
5  *
6  * Extracts the slow control data (everything except the primary event ID 1)
7  * from a MIDAS data file.
8  */
9 
10 #include <stdio.h>
11 #include <stdlib.h>
12 
13 #ifndef __CINT__
14 //#include "/home/nam/fadc_software/midas/include/midas.h"
15 #include "midas.h"
16 #include "../analyzer/compress/mucap_compress.h"
17 #endif
18 
19 #include "TH1.h"
20 #include "TCanvas.h"
21 
22 #define PRINT_SAMPLES 1
23 
24 TCanvas *c = 0;
25 void process_file(char *input_filename, int addr, int fadc);
26 int process_event(void *pevent, int addr, int fadc);
27 
28 
29 bool get_input_event(FILE * fp, char *pevent, int event_number)
30 {
31  int event_size, size_read;
32  EVENT_HEADER *header;
33 
34  size_read = fread(pevent, 1, sizeof(EVENT_HEADER), fp);
35  if (size_read != sizeof(EVENT_HEADER))
36  return !SUCCESS;
37 
38  header = (EVENT_HEADER *) pevent;
39  printf("block# %d\n", header->serial_number);
40 
41  if (event_number == 0 && header->trigger_mask != MIDAS_MAGIC) {
42  fprintf(stderr, "Input is not a valid MIDAS file\n");
43  exit(1);
44  }
45 
46  event_size = header->data_size;
47 
48  if ((unsigned int) event_size >= MAX_EVENT_SIZE - sizeof(EVENT_HEADER)) {
49  fprintf(stderr, "Invalid event size %d\n", event_size);
50  exit(1);
51  }
52 
53  size_read = fread(pevent + sizeof(EVENT_HEADER), 1, event_size, fp);
54  if (size_read != event_size)
55  return !SUCCESS;
56 
57  return SUCCESS;
58 }
59 
60 void process_file(char *input_filename, int addr, int fadc)
61 {
62  FILE *input_fp;
63  char *input_event;
64  int event_number;
65 
66  c = new TCanvas();
67 
68  char host_name[256], exp_name[32];
69  cm_get_environment(host_name, sizeof(host_name), exp_name, sizeof(exp_name));
70  cm_connect_experiment(host_name, exp_name, "show_pulses", NULL);
71  cm_get_experiment_database(&hDB, NULL);
73 
74  /* Open input and output files */
75  input_fp = fopen(input_filename, "rb");
76  if (!input_fp) {
77  fprintf(stderr, "Unable to open %s: ", input_filename);
78  perror("");
79  exit(1);
80  }
81 
82  /* Allocate buffer to hold MIDAS event */
83  input_event = (char *) malloc(MAX_EVENT_SIZE*2);
84 
85  /* Loop over MIDAS events */
86  event_number = 0;
87  while (get_input_event(input_fp, input_event, event_number) == SUCCESS) {
88 
89  EVENT_HEADER *header = (EVENT_HEADER *) input_event;
90 
91  if(header->event_id == 1) {
92  //expand_event(header+1, header+1);
93  int stop = process_event(header + 1, addr, fadc);
94  if(stop) break;
95  }
96 
97  event_number++;
98  printf("event_number %d\n", event_number);
99  }
100 
101  fclose(input_fp);
102 }
103 
104 int processIsland(int t0, int nsamples, int *samples, int channel)
105 {
106 
107  c->Clear();
108  char title[80];
109  sprintf(title, "Waveform: channel %d, t0=%d", channel, t0);
110 
111  double sum = 0;
112  int min = 4097;
113  int max = -1;
114 
115  double binWidth = 100./100.00;
116  TH1D *waveform = new TH1D("waveform", title, nsamples, -0.5*binWidth, (nsamples-0.5)*binWidth);
117  for(int i = 0; i < nsamples; i++) {
118  waveform->SetBinContent(i+1, samples[i]);
119  sum += samples[i];
120 
121  if(samples[i] < min) min = samples[i];
122  if(samples[i] > max) max = samples[i];
123  }
124  waveform->Draw();
125 
126  double mean = sum/nsamples;
127  double sumSqrDeviations = 0;
128 
129  for(int i = 0; i < nsamples; i++) {
130  double deviation = samples[i] - mean;
131  sumSqrDeviations += deviation*deviation;
132  }
133  double rms = sqrt(sumSqrDeviations/(nsamples-1));
134 
135  printf("Number of samples = %d, Min = %d, Max = %d, Mean = %lf, RMS = %lf\n", nsamples, min, max, mean, rms);
136 
137  // if(max > 3500) {
138  c->Update();
139  c->Modified();
140 
141  char dummy[80];
142  fgets(dummy, sizeof(dummy), stdin);
143  if(dummy[0] == 'q') {
144  return 1;
145  }
146 //}
147 
148  delete waveform;
149  //delete c;
150 
151 return 0;
152 }
153 
154 int process_event(void *pevent, int addr, int fadc)
155 {
156  for(int i = 0; i < 8; i++) {
157 
158  if(i != fadc && fadc >= 0) continue;
159 
160  // get the bank
161  char name[80];
162  sprintf(name, "N%c%02x", 'a' + i, addr);
163  unsigned char *raw;
164  int bank_size = bk_locate(pevent, name, &raw);
165  printf("bank_size %d,", bank_size);
166  bank_size = bank_size / 10;
167  printf(" npkt %d\n", bank_size);
168 
169 #ifdef PRINT_SAMPLES
170  printf("-------------------------------------------------------------------------------\n");
171 #endif
172 
173  int islandSamples[8192];
174  int islandNsamples = 0;
175  int islandT0 = -1;
176  bool firstIsland = true;
177 
178  // loop through words to build up "islands"
179  int lastTimestamp = 0;
180  for(int j = 0; j < bank_size; j++) {
181  printf("pkt %d\n", j);
182 
183  // data format:
184  //
185  // bits
186  // 78-52 timestamp
187  // 51-48 overflow
188  // 47-36 sampleB0
189  // 35-24 sampleA0
190  // 23-12 sampleB1
191  // 11-0 sampleA1
192  int timestamp = (raw[j*10+0] << 20) |
193  (raw[j*10+1] << 12) |
194  (raw[j*10+2] << 4) |
195  (raw[j*10+3] >> 4);
196  bool overflowB0 = ((raw[j*10+3] & 0x08) != 0);
197  bool overflowA0 = ((raw[j*10+3] & 0x04) != 0);
198  bool overflowB1 = ((raw[j*10+3] & 0x02) != 0);
199  bool overflowA1 = ((raw[j*10+3] & 0x01) != 0);
200  int sampleB0 = (overflowB0 << 12) |
201  (raw[j*10+4] << 4) |
202  (raw[j*10+5] >> 4);
203  int sampleA0 = (overflowA0 << 12) |
204  ((raw[j*10+5] & 0xf) << 8) |
205  (raw[j*10+6]);
206  int sampleB1 = (overflowB1 << 12) |
207  (raw[j*10+7] << 4) |
208  (raw[j*10+8] >> 4);
209  int sampleA1 = (overflowA1 << 12) |
210  ((raw[j*10+8] & 0xf) << 8) |
211  (raw[j*10+9]);
212 
213  if(timestamp < lastTimestamp) {
214  printf("Time ordering error!\n");
215  }
216 
217  if(timestamp != lastTimestamp + 1) {
218  if(!firstIsland) {
219  int stop = processIsland(islandT0*2, islandNsamples, islandSamples, i);
220 #ifdef PRINT_SAMPLES
221  printf("*\n");
222 #endif
223 
224  if(stop) {
225  return 1;
226  }
227  }
228  islandT0 = timestamp;
229  islandNsamples = 0;
230  firstIsland = false;
231  }
232  lastTimestamp = timestamp;
233 
234 #ifdef PRINT_SAMPLES
235  printf("fadc=%d t=%d %d %d %d %d\n", i, timestamp,
236  sampleA1, sampleB1, sampleA0,sampleB0);
237 #endif
238 
239  islandSamples[islandNsamples++] = sampleA1;
240  islandSamples[islandNsamples++] = sampleB1;
241  islandSamples[islandNsamples++] = sampleA0;
242  islandSamples[islandNsamples++] = sampleB0;
243  }
244 
245  if(islandNsamples != 0) {
246  int stop = processIsland(islandT0*2, islandNsamples, islandSamples, i);
247 #ifdef PRINT_SAMPLES
248  printf("*\n");
249 #endif
250  if (stop) return 1;
251  }
252  }
253 
254  return 0;
255 }