CoMPASS_BinReader/Bin2Root.cpp

#include "BinReader.h"

#include <sys/time.h> /** struct timeval, select() */

inline unsigned int getTime_us(){
  unsigned int time_us;
  struct timeval t1;
  struct timezone tz;
  gettimeofday(&t1, &tz);
  time_us = (t1.tv_sec) * 1000 * 1000 + t1.tv_usec;
  return time_us;
}

#include "TFile.h"
#include "TTree.h"
#include "TMacro.h"

#define MAX_TRACE_LENGTH 2000
#define MAX_MULTI  100
#define NMINARG 5
#define debug 0 // for > 1 number of hit to debug

#include <sstream>
#include <string>
#include <vector>
#include <utility>

std::vector<std::pair<int, int64_t>> readFileAndExtractData(const std::string& filename) {
  std::vector<std::pair<int, int64_t>> data;
  std::ifstream file(filename);
  
  if (!file) {
    std::cerr << "Error opening file!" << std::endl;
    return data;  // Return empty vector on failure
  }

  std::string line;
  while (std::getline(file, line)) {
    std::istringstream iss(line);
    int first;
    int64_t second;

    // Extract an int and an int64_t from the line
    if (iss >> first >> second) {
        data.emplace_back(first, second);
    } else {
        std::cerr << "Error parsing line: " << line << std::endl;
    }
  }

  file.close();
  return data;
}

int extractDigiSN(const std::string& str) {
    std::size_t firstUnderscore = str.find('_');
    if (firstUnderscore == std::string::npos) {
        return -1;  // No first underscore found, return an error value
    }

    std::size_t secondUnderscore = str.find('_', firstUnderscore + 1);
    if (secondUnderscore == std::string::npos) {
        return -1;  // No second underscore found, return an error value
    }

    // Extract the substring between the first and second underscores
    std::string extracted = str.substr(firstUnderscore + 1, secondUnderscore - firstUnderscore - 1);
    
    // Convert the extracted string to an integer using std::stoi
    try {
        return std::stoi(extracted);
    } catch (const std::invalid_argument& e) {
        std::cerr << "Invalid number format: " << extracted << std::endl;
        return -1;  // Return an error value on conversion failure
    } catch (const std::out_of_range& e) {
        std::cerr << "Number out of range: " << extracted << std::endl;
        return -1;  // Return an error value if the number is too large
    }
}

//^#############################################################
//^#############################################################
int main(int argc, char **argv) {
  
  printf("=========================================\n");
  printf("===      CoMPASS Binary to root       ===\n");
  printf("=========================================\n");  
  if (argc < NMINARG)    {
    printf("Incorrect number of arguments:\n");
    printf("%s [timeWindow] [withTrace] [time Offset File] [inFile1]  [inFile2] .... \n", argv[0]);
    printf("       timeWindow : in ns, -1 = no event building \n");   
    printf("        withTrace : 0 for no trace, 1 for trace \n");   
    printf(" time Offset File : 0 for nothing \n");   
    printf("\n");
    printf(" Example: %s -1 0  '\\ls -1 *001*.bin' (no event build, no trace, no verbose)\n", argv[0]);
    printf("          %s 100 0  '\\ls -1 *001*.bin' (event build with 100 ns, no trace, no verbose)\n", argv[0]);
    printf("\n");
    printf(" Time offset file format :\n");
    printf("[digi ID]   [time offset(int64_t)]\n");
    printf("\n");
    printf(" file name should be this format CH0@V1725_324_Data_run_196.bin\n");
    printf("     so that the serial number is between the 1st and 2nd underscore\n");
    printf("\n\n");

    return 1;
  }

  unsigned int runStartTime = getTime_us();

  ///============= read input
  long timeWindow = atoi(argv[1]);
  bool traceOn = atoi(argv[2]);
  std::string timeOffsetFile = argv[3];
  const int nFile = argc - NMINARG + 1;
  TString inFileName[nFile];
  for( int i = 0 ; i < nFile ; i++){ inFileName[i] = argv[i + NMINARG -1];}
  
  /// Form outFileName;
  int pos = inFileName[0].Last('_');
  TString outFileName = "";
  if( pos == -1 ) {
    if( timeWindow < 0 ) {
      outFileName ="temp_single.root"; 
    }else{
      outFileName = "temp_" + std::to_string(timeWindow) + ".root"; 
    }
  }else{
    TString temp = inFileName[0];
    temp.Remove(0, pos+1);
    int pos = temp.Last('.');
    temp.Remove(pos);
    outFileName = Form("run%03d", atoi(temp));
    if( timeWindow < 0 ) {
      outFileName +="_single.root"; 
    }else{
      outFileName += "_" + std::to_string(timeWindow) + ".root"; 
    }
  }

  std::vector<std::pair<int, int64_t>> timeOffsetList ;

  printf("-------> Out file name : %s \n", outFileName.Data());
  printf("=========================================\n"); 
  printf("       Time Window = %ld ns = %.1f us\n", timeWindow, timeWindow/1000.);
  printf("     Include Trace = %s\n", traceOn ? "Yes" : "No");
  printf("    Max multiplity = %d hits/event (hard coded)\n", MAX_MULTI);
  if( traceOn ) printf("  Max Trace Length = %d (hard coded)\n", MAX_TRACE_LENGTH);
  if( timeOffsetFile != "0" ) {
    printf("  Time Offset File = %s\n", timeOffsetFile.c_str());
    timeOffsetList = readFileAndExtractData(timeOffsetFile);
    for (size_t i = 0 ; timeOffsetList.size(); i++) {
      printf("Digi-%4d   %16llu ps\n", timeOffsetList[i].first, timeOffsetList[i].second);
    }
  }
  printf("========================================= Number of Files : %d \n", nFile);

  unsigned long long int totalHitCount = 0;

  BinReader reader[nFile];
  for(int i = 0; i < nFile; i++ ){
    printf(">>>>>>>> %2d | %s \n", i, inFileName[i].Data());
    reader[i].OpenFile(inFileName[i]);
    for (size_t i = 0 ; timeOffsetList.size(); i++) {
      int sn = extractDigiSN(inFileName[i].Data());
      if( sn == timeOffsetList[i].first ) reader[i].SetTimeOffset(timeOffsetList[i].second);
    }
    reader[i].ScanNumHit();
    totalHitCount += reader[i].GetNumHit();
  }

  printf("======================= total Hit Count : %llu\n", totalHitCount);

  printf("========================================= Initializing tree...\n"); 

  unsigned long long                evID = 0;
  unsigned int                     multi = 0;
  unsigned short           sn[MAX_MULTI] = {0}; 
  unsigned short           ch[MAX_MULTI] = {0}; 
  unsigned short            e[MAX_MULTI] = {0};
  unsigned short           e2[MAX_MULTI] = {0};
  unsigned long long      e_t[MAX_MULTI] = {0};
  unsigned short          e_f[MAX_MULTI] = {0};
  unsigned short  traceLength[MAX_MULTI];
  short trace[MAX_MULTI][MAX_TRACE_LENGTH];

  TFile * outRootFile = new TFile(outFileName, "recreate");
  TTree * tree = new TTree("tree", outFileName);
  tree->Branch("evID",           &evID, "event_ID/l"); 
  tree->Branch("multi",         &multi, "multi/i"); 
  tree->Branch("sn",                sn, "sn[multi]/s");
  tree->Branch("ch",                ch, "ch[multi]/s");
  tree->Branch("e",                  e, "e[multi]/s");
  tree->Branch("e2",                e2, "e2[multi]/s");
  tree->Branch("e_t",              e_t, "e_t[multi]/l");
  tree->Branch("e_f",              e_f, "e_f[multi]/s");
  tree->Branch("traceLength", traceLength, "traceLength[multi]/s");
  
  if( traceOn ) {
    tree->Branch("trace", trace,"trace[multi][MAX_TRACE_LENGTH]/S");
    tree->GetBranch("trace")->SetCompressionSettings(205);
  }

  printf("========================================= Event Building ...\n"); 

  //set TimeWindow to ps;
  timeWindow *= 1000;

  // record the first time stamp and last time stamp
  unsigned long long tStart = 0;
  unsigned long long tEnd = 0;

  unsigned long long t0 = -1;
  short g0 = 0 ;
  int nFileFinished = 0;
  multi = 0;
  evID = 0;

  std::vector<Data> events;
  long long ID[nFile]; // filled hit ID
  for( int i = 0 ; i < nFile ; i++ ) {
    reader[i].ReadBlock();
    ID[i] = -1;
  }

  unsigned long long hitProcessed = 0;

  do{

    if( debug ) printf("################################ ev build %llu \n", evID);
    events.clear();
    
    // //find earliest time group;
    t0 = -1;
    for( short i = 0; i < nFile; i++){
      if( debug ) printf("NNNNNNNNNNNNNNN %d | ID: %lld, %lld, %lld \n", i, ID[i], reader[i].GetHitID(), reader[i].GetNumHit() );
      if( ID[i] + 1 >= reader[i].GetNumHit())  continue;
      // if( ID[i] >= reader[i].GetHitID() ) reader[i].ReadBlock();
      if( debug ) reader[i].data.Print();
      if( reader[i].data.TimeStamp < t0 ) {
        t0 = reader[i].data.TimeStamp;
        g0 = i;
      }
    }

    if( evID == 0 ) tStart = t0;
    if( debug ) printf("First timestamp is %llu, file ID : %u\n", t0, g0);

    for( short i = 0; i < nFile; i++){

      if( ID[i] + 1 >= reader[i].GetNumHit() ) continue;
      
      if( timeWindow >= 0 ){
        do{
          if( (long int)( reader[i].data.TimeStamp - t0) <= timeWindow ){
            events.push_back(reader[i].data);
            ID[i] ++;
          }else{
            break;
          }
        }while( reader[i].ReadBlock() == 1  );
      }else{
        events.push_back(reader[g0].data);
        ID[g0] ++;
        reader[g0].ReadBlock();
        break;
      }

      //if( timeWindow < 0) break;
    }

    if( events.size() > 1 ){
      std::sort(events.begin(), events.end(), [](const Data& a, const Data& b) {
        return a.TimeStamp < b.TimeStamp;
      });
    }
    
    tEnd = events.back().TimeStamp;

    hitProcessed += events.size();    
    if( hitProcessed % (traceOn ? 10000 : 10000) == 0 ) printf("hit Porcessed %llu/%llu hit....%.2f%%\n\033[A\r", hitProcessed, totalHitCount, hitProcessed*100./totalHitCount);

    multi = events.size() ;
    if( events.size() >= MAX_MULTI ) {
      printf("\033[31m event %lld has size = %d > MAX_MULTI = %d \033[0m\n", evID, multi, MAX_MULTI);
      multi = MAX_MULTI;
    }
    if( debug ) printf("=================================== filling data | %u \n", multi);

    for( size_t p = 0; p < multi ; p ++ ) {
      if( debug ) {printf("%4zu | ", p); events[p].Print();}

      sn[p] = events[p].BoardID;
      ch[p] = events[p].Channel;
      e[p] = events[p].Energy;
      e2[p] = events[p].Energy_short;
      e_t[p] = events[p].TimeStamp / 1000;
      e_f[p] = events[p].TimeStamp % 1000;
      
      traceLength[p] = events[p].NSample;
      if( traceOn ){
        if( traceLength[p] > MAX_TRACE_LENGTH ) {
          printf("\033[31m event %lld has trace length = %d > MAX_TRACE_LENGTH = %d \033[0m\n", evID, traceLength[p], MAX_TRACE_LENGTH);
          traceLength[p] = MAX_TRACE_LENGTH;
        }

        for( int hh = 0; hh < traceLength[p]; hh++){
          trace[p][hh] = events[p].Trace[hh];
        }
      }
    }

    outRootFile->cd();
    tree->Fill();
    // tree->Write();

    multi = 0;
    evID ++;

    nFileFinished = 0;
    for( int i = 0; i < nFile; i++ ){
      // printf(" %d | %lld  %lld \n", i, ID[i], reader[i].GetNumHit() );
      if( ID[i] + 1 >= reader[i].GetNumHit() ) nFileFinished ++ ; 
    }

    if( debug ) printf(" >>>> nFileFinished : %d \n", nFileFinished);

    if( events.size() == 0 ) break;

    if( debug > 1 && hitProcessed > debug ) break;

  }while( nFileFinished < nFile );

  tree->Write();

  unsigned int runEndTime = getTime_us();
  double runTime = (runEndTime - runStartTime) * 1e-6;
  printf("######################################### finished.\n");
  printf(" event building time = %.2f sec = %.2f min\n", runTime, runTime/60.);
  // printf("  total events built = %llu by event builder (%llu in tree)\n", evID, tree->GetEntriesFast());
  printf(" total hit processed = %llu\n", hitProcessed);
  printf("  total events built = %llu by event builder\n", evID);
  double tDuration_sec = (tEnd - tStart) * 1e-12;
  printf("     first timestamp = %20llu ps\n", tStart);
  printf("      last timestamp = %20llu ps\n", tEnd);
  printf(" total data duration = %.2f sec = %.2f min\n", tDuration_sec, tDuration_sec/60.);
  printf("========================================> saved to %s \n", outFileName.Data());

  TMacro info;
  info.AddLine(Form("tStart= %20llu ns",tStart));
  info.AddLine(Form("  tEnd= %20llu ns",tEnd));
  info.Write("info");
  outRootFile->Close();

  return 0;

}