Timepix3EventLoader.C 16.9 KB
Newer Older
1
2
3
4
5
6
7
#include "Timepix3EventLoader.h"

#include <dirent.h>
#include <sstream>
#include <fstream>
#include <string>
#include <stdio.h>
Daniel Hynds's avatar
Daniel Hynds committed
8
#include <stdint.h>
9
#include <stdlib.h>
Daniel Hynds's avatar
Daniel Hynds committed
10
#include <bitset>
11
12
13
14
15

Timepix3EventLoader::Timepix3EventLoader(bool debugging)
: Algorithm("Timepix3EventLoader"){
  debug = debugging;
  applyTimingCut = false;
16
  m_currentTime = 0.;
17
  m_minNumberOfPlanes = 1;
18
  m_prevTime = 0;
19
  m_shutterOpen = false;
20
21
22
23
24
25
}


void Timepix3EventLoader::initialise(Parameters* par){
 
  parameters = par;
26
 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
  // Take input directory from global parameters
  m_inputDirectory = parameters->inputDirectory;
  
  // File structure is RunX/ChipID/files.dat
  
  // Open the root directory
  DIR* directory = opendir(m_inputDirectory.c_str());
  if (directory == NULL){tcout<<"Directory "<<m_inputDirectory<<" does not exist"<<endl; return;}
  dirent* entry; dirent* file;

  // Read the entries in the folder
  while (entry = readdir(directory)){
    
    // If these are folders then the name is the chip ID
    if (entry->d_type == DT_DIR){
      
      // Open the folder for this device
      string detectorID = entry->d_name;
      string dataDirName = m_inputDirectory+"/"+entry->d_name;
      DIR* dataDir = opendir(dataDirName.c_str());
      string trimdacfile;
      
49
50
      // Check if this device is to be masked
      if(parameters->masked.count(detectorID) != 0) continue;
51
52
53
54

      // Check if this device has conditions loaded and is a Timepix3
      if(parameters->detector.count(detectorID) == 0) continue;
      if(parameters->detector[detectorID]->m_detectorType != "Timepix3") continue;
55
      
56
57
58
59
60
      // Get all of the files for this chip
      while (file = readdir(dataDir)){
        string filename = dataDirName+"/"+file->d_name;

        // Check if file has extension .dat
CLICdp user's avatar
CLICdp user committed
61
        if (string(file->d_name).find("-1.dat") != string::npos){
62
63
64
65
66
67
          m_datafiles[detectorID].push_back(filename.c_str());
          m_nFiles[detectorID]++;

          // Initialise null values for later
          m_currentFile[detectorID] = NULL;
          m_fileNumber[detectorID] = 0;
68
69
          m_syncTime[detectorID] = 0;
          m_clearedHeader[detectorID] = false;
70
71
72
73
        }
        
        // If not a data file, it might be a trimdac file, with the list of masked pixels etc.
        if (string(file->d_name).find("trimdac") != string::npos){
74
75
          // If we have already found the masked trimdac file, use it for preference
          if(trimdacfile.find("masked") != string::npos) continue;
76
77
78
79
80
          trimdacfile = filename;
        }
        
      }
      
81
      // If files were stored, register the detector (check that it has alignment data)
82
      if(m_nFiles.count(detectorID) > 0){
83
84
        
        tcout<<"Registering detector "<<detectorID<<endl;
85
        if(parameters->detector.count(detectorID) == 0){tcout<<"Detector "<<detectorID<<" has no alignment/conditions loaded. Please check that it is in the alignment file"<<endl; return;}
86
//        parameters->registerDetector(detectorID);
87
88
89
90
91
92
93
94
95
96
97
98
99
        
        // Now that we have all of the data files and mask files for this detector, pass the mask file to parameters
        tcout<<"Set mask file "<<trimdacfile<<endl;
        parameters->detector[detectorID]->setMaskFile(trimdacfile);
        
        // Apply the pixel masking
        maskPixels(detectorID,trimdacfile);
      }

    }
  }
}

Daniel Hynds's avatar
Daniel Hynds committed
100
StatusCode Timepix3EventLoader::run(Clipboard* clipboard){
101
  
102
103
104
105
  // This will loop through each timepix3 registered, and load data from each of them. This can
  // be done in one of two ways: by taking all data in the time interval (t,t+delta), or by
  // loading a fixed number of pixels (ie. 2000 at a time)
  
106
//  tcout<<"== New event"<<endl;
107
  int endOfFiles = 0; int devices = 0; int loadedData = 0;
108

109
110
  // Loop through all registered detectors
  for(int det = 0; det<parameters->nDetectors; det++){
111
    
112
113
114
    // Check if they are a Timepix3
    string detectorID = parameters->detectors[det];
    if(parameters->detector[detectorID]->type() != "Timepix3") continue;
115
    
116
    // Make a new container for the data
117
    Pixels* deviceData = new Pixels();
118
119
    SpidrSignals* spidrData = new SpidrSignals();

120
		// Load the next chunk of data
121
122
    if(debug) tcout<<"Loading data from "<<detectorID<<endl;
    bool data = loadData(detectorID,deviceData,spidrData);
123
    
124
125
    // If data was loaded then put it on the clipboard
    if(data){
126
      loadedData++;
127
      if(debug) tcout<<"Loaded "<<deviceData->size()<<" pixels for device "<<detectorID<<endl;
128
129
      clipboard->put(detectorID,"pixels",(TestBeamObjects*)deviceData);
    }
130
    clipboard->put(detectorID,"SpidrSignals",(TestBeamObjects*)spidrData);
131
132
133
134
    
    // Check if all devices have reached the end of file
    devices++;
    if(feof(m_currentFile[detectorID])) endOfFiles++;
135
136
  }
  
137
138
139
  // Increment the event time
  parameters->currentTime += parameters->eventLength;
  
140
  // If all files are finished, tell the event loop to stop
Daniel Hynds's avatar
Daniel Hynds committed
141
  if(endOfFiles == devices) return Failure;
142
  
143
  // If no/not enough data in this event then tell the event loop to directly skip to the next event
Daniel Hynds's avatar
Daniel Hynds committed
144
  if(loadedData < m_minNumberOfPlanes) return NoData;
145
146
  
  // Otherwise tell event loop to keep running
147
  cout<<"\rCurrent time: "<<std::setprecision(4)<<std::fixed<<parameters->currentTime<<flush;
148
//  cout<<endl;
Daniel Hynds's avatar
Daniel Hynds committed
149
  return Success;
150
151
152
153
154
155
156
}

// Function to load the pixel mask file
void Timepix3EventLoader::maskPixels(string detectorID, string trimdacfile){
  
  // Open the mask file
  ifstream trimdacs;
Daniel Hynds's avatar
Daniel Hynds committed
157
  trimdacs.open(trimdacfile.c_str());
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
  
  // Ignore the file header
  string line;
  getline(trimdacs,line);
  int t_col, t_row, t_trim, t_mask, t_tpen;

  // Loop through the pixels in the file and apply the mask
  for(int col=0;col<256;col++){
    for(int row=0;row<256;row++){
      trimdacs >> t_col >> t_row >> t_trim >> t_mask >> t_tpen;
      if(t_mask) parameters->detector[detectorID]->maskChannel(t_col,t_row);
    }
  }
  
  // Close the files when finished
  trimdacs.close();

}

// Function to load data for a given device, into the relevant container
178
bool Timepix3EventLoader::loadData(string detectorID, Pixels* devicedata, SpidrSignals* spidrData){
179

180
181
182
//  if(detectorID == "W0019_F07") debug = true;
//  if(detectorID != "W0019_F07") debug = false;
  
183
  bool extra = false; //temp
184
//  if(detectorID == parameters->DUT) extra = true;
185
186
//  if(detectorID == "W0002_J05") extra = true;
  
187
  if(debug) tcout<<"Loading data for device "<<detectorID<<endl;
188

189
  // Check if current file is open
190
  if(m_currentFile[detectorID] == NULL || feof(m_currentFile[detectorID])){
191
  if(debug) tcout<<"No current file open "<<endl;
192
193
194

    // If all files are finished, return
    if(m_fileNumber[detectorID] == m_datafiles[detectorID].size()){
195
      if(debug) tcout<<"All files have been analysed. There were "<<m_datafiles[detectorID].size()<<endl;
196
197
198
      return false;
    }
    
199
200
    // Open a new file
    m_currentFile[detectorID] = fopen(m_datafiles[detectorID][m_fileNumber[detectorID]].c_str(),"rb");
201
    tcout<<"Loading file "<<m_datafiles[detectorID][m_fileNumber[detectorID]]<<endl;
202
   
203
204
    // Mark that this file is done
    m_fileNumber[detectorID]++;
205
206
    
    // Skip the header - first read how big it is
207
208
    uint32_t headerID;
    if (fread(&headerID, sizeof(headerID), 1, m_currentFile[detectorID]) == 0) {
209
      tcout << "Cannot read header ID for device " << detectorID << endl;
210
211
      return false;
    }
212
    
213
214
215
    // Skip the rest of the file header
    uint32_t headerSize;
    if (fread(&headerSize, sizeof(headerSize), 1, m_currentFile[detectorID]) == 0) {
216
      tcout << "Cannot read header size for device " << detectorID << endl;
217
218
      return false;
    }
219
    
220
221
222
223
224
225
226
227
		// Finally skip the header
    rewind(m_currentFile[detectorID]);
    fseek(m_currentFile[detectorID], headerSize, SEEK_SET);
  }
  
  // Now read the data packets.
  ULong64_t pixdata = 0; UShort_t thr = 0;
  int npixels=0;
228
229
  bool fileNotFinished = false;
  
230
231
  // Read till the end of file (or till break)
  while (!feof(m_currentFile[detectorID])) {
232
233
    
    // Read one 64-bit chunk of data
234
    const int retval = fread(&pixdata, sizeof(ULong64_t), 1, m_currentFile[detectorID]);
235
236
237
    if(debug) bitset<64> packetContent(pixdata);
    if(debug) tcout<<hex<<pixdata<<dec<<endl;
    if(debug) tcout<<pixdata<<endl;
238
    if (retval == 0) continue;
239
240
241
    
    // Get the header (first 4 bits) and do things depending on what it is
    // 0x4 is the "heartbeat" signal, 0xA and 0xB are pixel data
242
243
    const UChar_t header = ((pixdata & 0xF000000000000000) >> 60) & 0xF;
    
244
245
246
247
248
249
250
251
    // Use header 0x4 to get the long timestamps (called syncTime here)
    if(header == 0x4){
      
      // The 0x4 header tells us that it is part of the timestamp
      // There is a second 4-bit header that says if it is the most
      // or least significant part of the timestamp
      const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF;

252
253
      // This is a bug fix. There appear to be errant packets with garbage data - source to be tracked down.
      // Between the data and the header the intervening bits should all be 0, check if this is the case
254
255
      const UChar_t intermediateBits = ((pixdata & 0x00FF000000000000) >> 48) & 0xFF;
      if(intermediateBits != 0x00) continue;
256
      
257
      // 0x4 is the least significant part of the timestamp
258
259
260
261
      if(header2 == 0x4){
        // The data is shifted 16 bits to the right, then 12 to the left in order to match the timestamp format (net 4 right)
        m_syncTime[detectorID] = (m_syncTime[detectorID] & 0xFFFFF00000000000) + ((pixdata & 0x0000FFFFFFFF0000) >> 4);
      }
262
			// 0x5 is the most significant part of the timestamp
263
264
265
266
267
268
      if(header2 == 0x5){
        // The data is shifted 16 bits to the right, then 44 to the left in order to match the timestamp format (net 28 left)
        m_syncTime[detectorID] = (m_syncTime[detectorID] & 0x00000FFFFFFFFFFF) + ((pixdata & 0x00000000FFFF0000) << 28);
        if( m_syncTime[detectorID] < 0x0000010000000000 && !m_clearedHeader[detectorID]) m_clearedHeader[detectorID] = true;
      }
      
269
//      if(extra) tcout<<"Updating heartbeat. Now syncTime = "<<(double)m_syncTime[detectorID]/(4096. * 40000000.)<<endl;
270
271
    }
    
272
273
274
    // Header 0x06 and 0x07 are the start and stop signals for power pulsing
    if(header == 0x0){
      
275
276
277
      // Only want to read these packets from the DUT
      if(detectorID != parameters->DUT) continue;

278
279
280
      // Get the second part of the header
      const UChar_t header2 = ((pixdata & 0x0F00000000000000) >> 56) & 0xF;

281
      // New implementation of power pulsing signals from Adrian
282
      if(header2 == 0x6){
283
        const uint64_t time( (pixdata & 0x0000000FFFFFFFFF) << 12 );
284
        
285
        const uint64_t controlbits = ((pixdata & 0x00F0000000000000) >> 52) & 0xF;
286

287
288
        const uint64_t powerOn = ((controlbits & 0x2) >> 1);
        const uint64_t shutterClosed = ((controlbits & 0x1));
289
290
        
        
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
        // Ignore packets if they arrive before the current event window
        if( parameters->eventLength != 0. && ((double)time/(4096. * 40000000.)) < (parameters->currentTime) ){
          continue;
        }
        
        // Stop looking at data if the signal is after the current event window (and rewind the file
        // reader so that we start with this signal next event)
        if( parameters->eventLength != 0. &&
           ((double)time/(4096. * 40000000.)) > (parameters->currentTime + parameters->eventLength) ){
          fseek(m_currentFile[detectorID], -1 * sizeof(ULong64_t), SEEK_CUR);
          fileNotFinished = true;
          break;
        }

        
        
        
        SpidrSignal* powerSignal = (powerOn ? new SpidrSignal("powerOn",time) : new SpidrSignal("powerOff",time));
        spidrData->push_back(powerSignal);
        if(debug) tcout<<"Power is "<< (powerOn ? "on" : "off") <<" power! Time: "<<std::setprecision(10)<<(double)time/(4096. * 40000000.)<<endl;
        
//				cout<<endl;
//				tcout<<"Shutter closed: "<<hex<<shutterClosed<<dec<<endl;
        
        SpidrSignal* shutterSignal = (shutterClosed ? new SpidrSignal("shutterClosed",time) : new SpidrSignal("shutterOpen",time));
        if(!shutterClosed){
          spidrData->push_back(shutterSignal);
          m_shutterOpen = true;
//          tcout<<"Have opened shutter with signal "<<shutterSignal->type()<<" at time "<<(double)time/(4096. * 40000000.)<<endl;
        }
        
        if(shutterClosed && m_shutterOpen){
          spidrData->push_back(shutterSignal);
          m_shutterOpen = false;
//          tcout<<"Have closed shutter with signal "<<shutterSignal->type()<<" at time "<<(double)time/(4096. * 40000000.)<<endl;
        }
327
        
328
        if(debug) tcout<<"Shutter is "<< (shutterClosed ? "closed" : "open") <<". Time: "<<std::setprecision(10)<<(double)time/(4096. * 40000000.)<<endl;
329
        
330
       }
331
332
333
334
335
336
337
338
339
      
      /*
      // 0x6 is power on
      if(header2 == 0x6){
        const uint64_t time( (pixdata & 0x0000000FFFFFFFFF) << 12 );
        SpidrSignal* signal = new SpidrSignal("powerOn",time);
        spidrData->push_back(signal);
        if(debug) tcout<<"Turned on power! Time: "<<(double)time/(4096. * 40000000.)<<endl;
      }
340
341
      // 0x7 is power off
      if(header2 == 0x7){
342
343
344
345
        const uint64_t time( (pixdata & 0x0000000FFFFFFFFF) << 12 );
        SpidrSignal* signal = new SpidrSignal("powerOff",time);
        spidrData->push_back(signal);
        if(debug) tcout<<"Turned off power! Time: "<<(double)time/(4096. * 40000000.)<<endl;
346
      }
347
       */
348
349
350
351
352
    }
    
    // In data taking during 2015 there was sometimes still data left in the buffers at the start of
    // a run. For that reason we keep skipping data until this "header" data has been cleared, when
    // the heart beat signal starts from a low number (~few seconds max)
353
354
    if(!m_clearedHeader[detectorID]) continue;
    
355
    // Header 0xA and 0xB indicate pixel data
356
    if(header == 0xA || header == 0xB) {
357
358
      
      // Decode the pixel information from the relevant bits
359
360
361
362
363
      const UShort_t dcol = ((pixdata & 0x0FE0000000000000) >> 52);
      const UShort_t spix = ((pixdata & 0x001F800000000000) >> 45);
      const UShort_t pix  = ((pixdata & 0x0000700000000000) >> 44);
      const UShort_t col = (dcol + pix / 4);
      const UShort_t row = (spix + (pix & 0x3));
364
      
365
366
367
      // Check if this pixel is masked
      if(parameters->detector[detectorID]->masked(col,row)) continue;

368
369
370
371
      // Get the rest of the data from the pixel
      const UShort_t pixno = col * 256 + row;
      const UInt_t data = ((pixdata & 0x00000FFFFFFF0000) >> 16);
      const unsigned int tot = (data & 0x00003FF0) >> 4;
372
373
374
      const uint64_t spidrTime(pixdata & 0x000000000000FFFF);
      const uint64_t ftoa(data & 0x0000000F);
      const uint64_t toa((data & 0x0FFFC000) >> 14);
375
      
376
      // Calculate the timestamp.
377
      long long int time = (((spidrTime << 18) + (toa << 4) + (15 - ftoa)) << 8) + (m_syncTime[detectorID] & 0xFFFFFC0000000000);
378
379
      if(debug) tcout<<"Pixel time "<<(double)time/(4096. * 40000000.)<<endl;
      if(debug) tcout<<"Sync time "<<(double)m_syncTime[detectorID]/(4096. * 40000000.)<<endl;
380
      
381
      // Add the timing offset from the coniditions file (if any)
382
      time += (long long int)(parameters->detector[detectorID]->timingOffset() * 4096. * 40000000.);
383

384
385
386
387
      // The time from the pixels has a maximum value of ~26 seconds. We compare the pixel time
      // to the "heartbeat" signal (which has an overflow of ~4 years) and check if the pixel
      // time has wrapped back around to 0
      
388
      // If the counter overflow happens before reading the new heartbeat
389
390
391
392
393
394
395
396
397
//      while( abs(m_syncTime[detectorID]-time) > 0x0000020000000000 ){
      if(!extra){
        while( m_syncTime[detectorID]-time > 0x0000020000000000 ){
          time += 0x0000040000000000;
        }
      }else{
        while( m_prevTime - time > 0x0000020000000000){
          time += 0x0000040000000000;
        }
398
399
      }

400
      // If events are loaded based on time intervals, take all hits where the time is within this window
401
402
403
      
      // Ignore pixels if they arrive before the current event window
      if( parameters->eventLength != 0. && ((double)time/(4096. * 40000000.)) < (parameters->currentTime) ){
404
405
406
        continue;
      }
      
407
408
      // Stop looking at data if the pixel is after the current event window (and rewind the file
      // reader so that we start with this pixel next event)
409
410
411
      if( parameters->eventLength != 0. &&
         ((double)time/(4096. * 40000000.)) > (parameters->currentTime + parameters->eventLength) ){
        fseek(m_currentFile[detectorID], -1 * sizeof(ULong64_t), SEEK_CUR);
412
        fileNotFinished = true;
413
414
415
        break;
      }
      
416
      // Otherwise create a new pixel object
417
      Pixel* pixel = new Pixel(detectorID,row,col,(int)tot,time);
418
419
      devicedata->push_back(pixel);
      npixels++;
420
      m_prevTime = time;
421
422

    }
423
    
424
    // Stop when we reach some large number of pixels (if events not based on time)
425
426
427
428
    if(parameters->eventLength == 0. && npixels == 2000){
      fileNotFinished = true;
      break;
    }
429
430
431
    
  }
  
432
  // If no data was loaded, return false
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
  if(npixels == 0) return false;
  
  return true;

}

void Timepix3EventLoader::finalise(){
  
}