/home/alex/Work/Research/Km3NeT/PPM-DU/ppm-du-v2.3/rawReader.cc

Go to the documentation of this file.

#include <iostream>
#include <fstream>
#include <string>
#include <iomanip>
#include <boost/lexical_cast.hpp>
#include <climits>
#include <vector>
#include <math.h>

#include "rawReader.h"
#include "domHeader.h"
#include "pmt.h"
#include "signal.h"
#include "constants.h"
//#include <TRandom3.h>

std::vector<std::vector<MyEvent> > allHits;

uint16_t
hit::Channel() {
  if (data == 0)
    throw "Hit un-initialized";
  uint8_t header = data[0];
  uint16_t chan = ((header >> 0) & 0x3); // 2 lower bits  of channel 
  chan |= ((header >> 3) & (0x7 << 2));  // 3 higher bits of channel
  return chan;
}

uint16_t
hit::Type() {
  if (data == 0)
    throw "Hit un-initialized"; 
  uint8_t header = reinterpret_cast<uint16_t*>(data)[0];
  uint16_t type = ((header >> 2) & 0x7); // bits 2-4
  return type;
}

uint32_t
hit::TimeStamp() {
  if (data == 0)
    throw "Hit un-initialized";
  uint32_t ts = reinterpret_cast<uint32_t*>(data)[0];
  utl::BinIO::swap(ts);
  return ts & 0x00FFFFFF;
}

uint16_t
hit::Pattern() {
  if (data == 0)
    throw "Hit un-initialized";
  uint16_t patt = reinterpret_cast<uint16_t*>(data + 4)[0];
  utl::BinIO::swap(patt);
  return patt;
}

void
hit::HexDump() {
  if (data == 0)
    throw "Hit un-initialized";
  for (size_t i = 0; i < SIZE_BYTES; ++i) {
    std::cout << std::hex << std::setfill('0') << std::setw(2) << (uint32_t) data[i];
  }
  std::cout << std::dec << std::endl;
}

void
hit::Dump() {
  if (data == 0)
    throw "Hit un-initialized"; 
  std::cout << std::dec
            << std::setfill(' ') << std::setw(3) << Channel()
            << std::setw(2) << Type()
            << std::setw(9) << TimeStamp() 
            << std::hex << std::setfill('0') << " 0x" << std::setw(4) << Pattern()
            << std::dec << std::endl;
}

/********************************************/

bool
frame::Read(std::istream& input) {
  if (rawPreamble == 0)
    rawPreamble = new uint8_t[PREAMBLE_SIZE_BYTES];
  input.read(reinterpret_cast<char*>(rawPreamble), PREAMBLE_SIZE_BYTES);
  if (input.fail())
    return false;
  if (rawData != 0)
    delete [] rawData;
  rawData = new uint8_t[DataSize_bytes()];
  input.read(reinterpret_cast<char*>(rawData), DataSize_bytes());
  if (input.fail())
    return false;
  return true;
}

uint8_t
frame::Status() {
  uint8_t st = GetNumber<uint8_t>(20u);
  return st;
}

uint32_t
frame::NbItems() {
  uint32_t Nb = reinterpret_cast<uint32_t*>(rawPreamble + 20u)[0];
  utl::BinIO::swap(Nb);
  return Nb & 0x00FFFFFF;
}

hit
frame::CurHit(size_t iHit) {
  if (rawData == 0)
    throw "frame not initialized";
  if (iHit >= NbItems())
    throw "Hit is out of bounds";
  return hit(rawData + iHit * hit::SIZE_BYTES);
}

void
frame::HexDumpPreamble() {
  if (rawPreamble == 0)
    throw "Raw Preamble un-initialized";
  std::cout << std::string(95, '=') << std::endl;
  for (size_t i = 0; i < PREAMBLE_SIZE_BYTES; ++i) {
    std::cout << std::hex << std::setfill('0') << std::setw(2) << (uint32_t) rawPreamble[i] << ' ';
  }
  std::cout << std::dec << std::endl;
}

void
frame::DumpPreamble() {
  if (rawPreamble == 0)
    throw "Raw Preamble un-initialized";
  std::cout << std::string(95, '=') << std::endl;
  std::cout << std::dec
            << std::setfill(' ') << std::setw(8) << FrameSize() 
            << std::setw(2) << ItemType()
            << std::setw(3) << FrameTarget()
            << std::setfill('0') << std::setw(17) << FrameTime()
            << std::setfill(' ') << std::setw(7) << FrameIndex()
            << std::setfill('0') << " 0x" << std::hex << std::setw(2) << Status() << std::dec
            << std::setfill(' ') << std::setw(9) << NbItems()
            << std::setw(6) << DomId()
            << std::setw(2) << StreamId()
            << std::setw(6) << RunNumber()
            << std::endl;
}

void
frame::HexDumpHits() {
  for (size_t i = 0; i < NbItems(); ++i) {
    std::cout << std::dec << std::setfill(' ') << std::setw(6) << i << ": ";
    CurHit(i).HexDump();
  }
  std::cout << std::dec << std::endl;
}

void
frame::DumpHits() {
  for (size_t i = 0; i < NbItems(); ++i) {
    std::cout << std::dec << std::setfill(' ') << std::setw(6) << i << ": ";
    CurHit(i).Dump();
  }
  std::cout << std::dec << std::endl;
}

void
frame::HexDumpTail() {
  if (rawData == 0)
    throw "Raw Data buffer un-initialized";
  for (size_t i = NbItems() * hit::SIZE_BYTES; i < DataSize_bytes(); ++i) {
    std::cout << std::hex << std::setfill('0') << std::setw(2) << (uint32_t) rawData[i] << ' ';
  }
  std::cout << std::dec << std::endl;
}

//**********************************************/

int cpt = 0;
sliceHead Header;

int
ReplacePmtId(const int dom, const int first) {
  if (dom == 103)
    return old2NewPmtMapCorrected2[first];
  return old2NewPmtMapCorrected01[first];
}

void
WriteSlice(dom& curDom) {
  //std::cout << " runNumber " << Header.runNumber << " slice " << Header.sliceIndex
  //        << " time " << Header.sliceTime << std::endl;
  domHeader curHeader;
  curHeader.SetRunNumber(Header.runNumber);
  curHeader.SetDomId(Header.domId);
  curHeader.SetDomIndex(Header.sliceIndex);
  curHeader.SetDomTime(Header.sliceTime);
  std::vector<pmt> pmtData;
  curDom.SetDomHeader(curHeader);
  //TRandom3 rand;  
  for (int nPm = 0; nPm < pmtNb; ++nPm) {
    if (!allHits.size()) continue;
    std::vector<MyEvent> event;
    event = allHits.at(nPm);
    int eventSize = event.size();
    //const int timeOffset = nPm*100;
    const int timeOffset = 0;
    pmt curPmt(ReplacePmtId(Header.domId, nPm), timeOffset);
    std::vector<signal> pmtSignals;
    if (!eventSize) continue;
    for (int i = 0; i < eventSize; ++i) {
      const long long int value = 8*event[i].timeStamp + event[i].firstBit + timeOffset;
      //const long long int value = rand.Uniform(pow(2, 27));
      //if ((nPm == 0) && (Header.domId == 101))
      //if (Header.domId == 101)
      //  std::cout << nPm << ' ' << 8*event[i].timeStamp + event[i].firstBit << ' ' << event[i].firstBit
      //            << ' ' << event[i].widthPat << std::endl;
      pmtSignals.push_back(signal(value, event[i].widthPat));
      //if (value < 0) std::cout << value << std::endl;
    }
    curPmt.SetSignal(pmtSignals);    
    pmtData.push_back(curPmt);
    pmtSignals.clear();
    curPmt.ClearSignalData();
    event.clear();
    event.resize(1);
    cpt += eventSize;
  }
  curDom.SetPmtsData(pmtData);
  pmtData.clear();      
}

int
Analyse(frame& thisFrame, bool toWrite) {
  int ichan;
  uint32_t curTS;
  int curFirstBit;
  int curLastBit;
  int curWidthPat;
  uint32_t lastTS[pmtNb];
  // event est un vecteur de MyEvent
  std::vector<MyEvent> event;
  MyEvent hitChan;
  uint16_t curPattern;
  int nbItems = thisFrame.NbItems();
  //int nbItems = int((thisFrame.FrameSize() - 9)*4./6.);
  //if (nbItems == 197379) return -1;
  //if (nbItems == 0) return -1;
  for(int i = 0; i < pmtNb; ++i) {
    lastTS[i] = 0; 
  } 
  allHits.resize(pmtNb);
  for (int it = 0; it < nbItems; ++it) {
    if(thisFrame.CurHit(it).data == 0) {
      return -1;
    }
    ichan = thisFrame.CurHit(it).Channel();
    curTS = thisFrame.CurHit(it).TimeStamp();
    curPattern = thisFrame.CurHit(it).Pattern();
    if((curFirstBit = GetFirstBit(curPattern, 0)) == -10)
      return -2;
    if((curLastBit = GetLastBit(curPattern, curFirstBit)) == -10)
      return -3;
    curWidthPat = curLastBit - curFirstBit;
    
    // c'est la que je bidouille
    if ((ichan > 30) || (ichan < 0)) continue;
    if (abs((int)(lastTS[ichan] - curTS)) == 2) {       
      if (allHits.at(ichan).size() > 0) {
        int previous = allHits.at(ichan).size() - 1;
        allHits.at(ichan).at(previous).widthPat += curWidthPat;
      }
    } else {
      hitChan.timeStamp = curTS;
      hitChan.firstBit = curFirstBit;
      hitChan.widthPat = curWidthPat;
      //if (ichan == 0)
      //  std::cout << ichan << ' ' << curTS << std::endl;
      // fill vector with structure
      allHits.at(ichan).push_back(hitChan);
    }
    // save time stamp and hit number of the channel
    lastTS[ichan] = curTS;      
    //if ((thisFrame.DomId() == 101) && (ichan == 1) && (nbItems < 100000))
  }
  /*
  if (toWrite == true) {
    if (thisFrame.DomId() == 102) {
      for (unsigned int ch = 0; ch < allHits.size(); ++ch) {
        for (unsigned int it = 0; it < allHits[ch].size(); ++it) {
          const uint64_t time = uint64_t((thisFrame.FrameIndex() - 1)*pow(2, 27)) + uint64_t(timeOffsetDom[int(thisFrame.DomId()) - 101]*1e9 + fineDomOffset);
          std::cout << time << ' ' << ch << ' '
                    << 8*allHits[ch][it].timeStamp + int(allHits[ch][it].firstBit) << ' '
                    << int(allHits[ch][it].widthPat) << std::endl;
        }
      }
    }
  }
  */
  return 1;
}

int
GetFirstBit(size_t pattern, int firstBit) {
  for (size_t i = firstBit; i < (int(sizeof(pattern))*CHAR_BIT); ++i) {
    if ((pattern & (1 << i)) != 0) {
      return i;
    }
  }
  return -10;
}

int
GetLastBit(size_t pattern, int firstBit) {
  for (int i = firstBit; i < (int(sizeof(pattern))*CHAR_BIT); ++i) {
    if ((pattern & (1 << i)) == 0) {
      return i;
    }
  }
  return -10;
}

long long int
FillRaw(std::ifstream& str, std::vector<dom>& vSlice, long long int bitPosition, bool toWrite) {
  cpt = 0;
  size_t nFrames = ~0u;
  frame curFrame;
  int nanalysed = 0;
  str.seekg(bitPosition);
  long long int newPosition = bitPosition;
  while (nFrames-- > 0) {
    bool toFill = true;
    if (not curFrame.Read(str))
      break;
    if (curFrame.StreamId() == 6) continue;
    Header.runNumber = curFrame.RunNumber();
    Header.domId = curFrame.DomId();
    Header.sliceIndex = curFrame.FrameIndex();
    const int fineDomOffset = fineTimeOffsetDom[curFrame.DomId() - 101];
    Header.sliceTime = uint64_t((curFrame.FrameIndex() - 1)*pow(2, 27) + timeOffsetDom[curFrame.DomId() - 101]*1e9) + fineDomOffset;
    //Header.sliceTime = uint64_t(curFrame.FrameTime());    
    if (toWrite == true) {
      std::cout << curFrame.FrameSize() << " 0 0 " << curFrame.FrameIndex() << ' '
              << (curFrame.FrameIndex() - 1 )*pow(2, 27) + timeOffsetDom[curFrame.DomId() - 101]*1e9 + fineDomOffset
              //<< 8*curFrame.FrameTime()
              << " 0 0 " << curFrame.FrameSize() << ' ' << curFrame.NbItems() << ' '
              << curFrame.DomId() << ' ' << curFrame.StreamId() << " 0" << std::endl;
    }
    
    allHits.clear();
    allHits.resize(0);
    switch(Analyse(curFrame, toWrite)) {
    case 0:
      toFill = false;
      break;
    case 1:
      toFill = true;
      break;
    case -1:
      std::cout << "error: empty frame" << std::endl;
      toFill = false;
      break;
    case -2:
      std::cout << "error: pb with firstBit" << std::endl;
      toFill = false;
      break;
    case -3:
      std::cout << "error: pb with lastBit" << std::endl;
      toFill = false;
      break;
    default:
      toFill = false;
      break;
    }
    /*
    if (toWrite == true) {
      if (curFrame.DomId() == 102) {
        for (int ch = 0; ch < allHits.size(); ++ch) {
          for (int it = 0; it < allHits[ch].size(); ++it) {
            std::cout << ch << ' ' << allHits[ch][it].timeStamp << ' '
                      << allHits[ch][it].firstBit << ' '
                      << allHits[ch][it].widthPat << std::endl;
          }
        }
      }
    }
    */
    nanalysed++;
    newPosition = str.tellg();
    if ((newPosition - bitPosition) > 2.0e9)
      break;
    if (toFill) {
      dom curDom;
      WriteSlice(curDom);
      //if (int(curFrame.NbItems()) == 197379) {
      //std::cout << int(curFrame.Status()) << std::endl;
      vSlice.push_back(curDom);
      //}
    }
  }
  std::cerr << "number of frames analysed: "<< nanalysed
            << " - total number of hits " << cpt << std::endl;
  return newPosition;
}

long long int 
GetFileLength(std::ifstream& str) {
  str.seekg(0, std::ifstream::end);
  long long int theSize = str.tellg();
  str.seekg(0, std::ifstream::beg);
  return theSize; 
}

bool
IsNoisySlice(const int signalNb, const double limitNb)
{
  return (signalNb > limitNb) ? true : false;
}