XIAEventBuilder/armory/EventBuilder.cpp

/*==================

Thie event builder both sort and build event. skip the *.to file.

===================*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cmath>
#include <stdbool.h>
#include <iostream>
#include <fstream>
#include <sys/time.h> /** struct timeval, select() */
#include <chrono>
#include <sstream>
#include <algorithm> // for std::remove_if
#include <cctype> // for std::isspace

#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TBenchmark.h"
#include "TStopwatch.h"
#include "TTreeIndex.h"

std::string trimSpaces(const std::string& str);
std::vector<std::string> split(const std::string& str, char delimiter);
unsigned int getTime_us();
unsigned long long getTime_ns();

#include "evtReader.h"

#define MAXMULTI 100
#define BUFFERSIZE 1000000 // number of time and filePos in buffer
#define DEBUG 0

int main(int argn, char **argv){
  printf("=====================================\n");
  printf("===   Event Builder from *.evt    ===\n");
  printf("=====================================\n");  
    
  if (argn < 6  )    {
    printf("Usage :\n");
    printf("%s [timeWindows] [mapping file] [Reject Flag] [QDC Flag] [SaveFileName] [*.evt File1] [*.evt File2] ...\n", argv[0]);
    printf("         timeWindows [int]: 1 unit = 10 ns \n");
    printf("   mapping file path [str]: the path of mapping file. \n");
    printf("         Reject Flag [int]: 0 = no rejection. see mapping.h\n");
    printf("                            1 = reject BGO\n");
    printf("                            2 = reject no gamma\n");
    printf("                            4 = reject no GAGG\n");
    printf("                            8 = reject zero energy data-point\n");
    printf("                            3 = reject BGO + no gamma, 5 = reject BGO + no GAGG, etc.\n");
    printf("            QDC Flag [int]: 0 = no qdc, 1 = with qdc\n");
    printf("        SaveFileName [str]: custom save file name \n");
    return 1;
  }

  int timeWindow = atoi(argv[1]);
  TString mappingFilePath = argv[2];
  unsigned short rejectFlag = atoi(argv[3]);
  unsigned short   qdcFlag = atoi(argv[4]);
  TString outFileName = argv[5];

  std::vector<std::string> inFileList;
  for( int i = 6; i < argn; i++) inFileList.push_back(argv[i]);

  printf(" Mapping file Path : %s \n", mappingFilePath.Data());
  printf("       Time window : %d ticks\n", timeWindow);
  printf("       Reject Flag : %u \n", rejectFlag);
  printf("      outfile Name : %s \n", outFileName.Data());

  printf("--------------- Number of in files %ld \n", inFileList.size());
  for( size_t i = 0; i < inFileList.size(); i++){
    printf("%2ld | %s \n", i, inFileList[i].c_str());
  }

  if( rejectFlag > 0 ) printf("================================== Rejection Filter Conditions\n");
  if( ( rejectFlag & 0x1 ) ) printf("\033[31m !!!! Reject event with BGO !!!! \033[0m\n");
  if( ( rejectFlag & 0x2 ) ) printf("\033[31m !!!! Reject event w/o Clover !!!! \033[0m\n");
  if( ( rejectFlag & 0x4 ) ) printf("\033[31m !!!! Reject event w/o GAGG !!!! \033[0m\n");
  if( ( rejectFlag & 0x8 ) ) printf("\033[31m !!!! Reject Zero-Energy Data Point !!!! \033[0m\n");
  if( rejectFlag > 0 ) printf(" Rejection filter does not apply to the last event when timeWindow >= 0 \n");
  printf("================================== Digesting Mapping file\n");

  std::ifstream mapFile( mappingFilePath.Data() );

  int mapping[16*12]; // fixed to be 12 digitizers
  for( int i = 0; i < 16*12 ; i++ ) mapping[i] = -1;

  if( !mapFile.is_open() ){
    printf("Cannot open mapping file : %s. Skip. \n", mappingFilePath.Data());
  }else{

    int index = 0;
    bool startMap = false;
    std::string line;
    while( std::getline(mapFile, line) ){
      // printf("|%s|\n", line.c_str());
      if( line.find("//-") != std::string::npos ) continue;
      if( line.find("<--") != std::string::npos ){
        startMap = true;
        continue;
      }

      if( startMap && line.find("<--") != std::string::npos ) break;
      if( startMap ){
        std::vector<std::string> list = split(line, ',');
        for( size_t k = 0; k < list.size(); k ++ ){
          if( list[k].find("//") != std::string::npos || k >= 16 ) continue;
          // printf("%ld | %s \n", k, list[k].c_str()); 
          mapping[index] = atoi(list[k].c_str());
          index ++;
        }

      }
    }
    mapFile.close();
  }

  //----  print mapping;
  for( int i = 0; i < 16*11; i++){
    if( i % 16 == 0 ) printf("Mod-%02d | ", i/16);
    if(                      mapping[i] <   0 ) printf("%4d,", mapping[i]);
    if(   0 <= mapping[i] && mapping[i] < 100 ) printf("\033[31m%4d\033[0m,", mapping[i]);
    if( 100 <= mapping[i] && mapping[i] < 200 ) printf("\033[34m%4d\033[0m,", mapping[i]);
    if( 200 <= mapping[i] && mapping[i] < 300 ) printf("\033[32m%4d\033[0m,", mapping[i]);
    if( 300 <= mapping[i] && mapping[i] < 400 ) printf("\033[35m%4d\033[0m,", mapping[i]);
    if( i % 16 == 15 ) printf("\n");
  }


  printf("================================== Creating Tree\n");

  printf(">>> Create output tree\n");
  TFile * saveFile = new TFile(outFileName, "recreate");
  TTree * newtree = new TTree("tree", outFileName);

  UInt_t eventID = 0 ;  
  UInt_t multi = 0; /// this is total multipicilty for all detectors
  newtree->Branch("multi", &multi, "multi/i");  
  newtree->Branch("evID", &eventID, "event_ID/i"); 
  
  // Int_t multiCry = 0 ; /// thi is total multiplicity for all crystal
  // newtree->Branch("multiCry", &multiCry, "multiplicity_crystal/I");  
  
  UInt_t     id[MAXMULTI] = {0};
  Int_t       e[MAXMULTI] = {-1};    
  ULong64_t e_t[MAXMULTI] = {0};
  UInt_t    qdc[MAXMULTI][8] = {0};
  newtree->Branch("id",       id, "id[multi]/i" );
  newtree->Branch("e",         e, "e[multi]/I" );
  newtree->Branch("e_t",     e_t, "e_timestamp[multi]/l");
  if( qdcFlag ) newtree->Branch("qdc",     qdc, "qdc[multi][8]/I");


  saveFile->cd();

  printf("================== Start processing....\n");
  Float_t Frac = 0.05; ///Progress bar
  TStopwatch StpWatch;
  StpWatch.Start();
  
  std::vector<timePos> hitList;
  evtReader * reader = nullptr;
  std::vector<DataBlock> event;
  event.clear();
  unsigned long int totalBlock;
  unsigned long int blockCount = 0 ;

  unsigned long long tStart = 0;
  unsigned long long tEnd = 0;

  unsigned int runStartTime = getTime_us();

  for( size_t i = 0 ; i < inFileList.size(); i++){

    reader = new evtReader(inFileList[i]);
    reader->ScanNumberOfBlock();
    totalBlock = reader->GetNumberOfBlock();
    blockCount = 0;

    do{

      hitList = reader->ReadBatchPos(BUFFERSIZE, DEBUG);
      blockCount += hitList.size();
      // std::cout << "please wait.....";
      printf("File-%ld %10lu block %10lu / %lu [%4.1f%%] | %u \r", i, hitList.size(), blockCount, totalBlock, blockCount*100./totalBlock, eventID);
      fflush(stdout);
      if( hitList.size() == 0 ) break;

      for( size_t k = 0; k < hitList.size(); k++ ){
        
        reader->ReadBlockAtPos(hitList[k].inFilePos);

        if( eventID == 0 ) tStart = reader->data->time;

        if( event.size() == 0 ) {
                    
          event.push_back(*(reader->data));

          if( timeWindow < 0 ) { // no event build
            multi = 1;
            int index  = event[0].crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (event[0].slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + event[0].ch;
            id[0] = mapping[index];
            e[0]   = event[0].energy;
            e_t[0] = event[0].time;

            if( qdcFlag ) for( int i = 0; i < 8; i++) qdc[0][i] = event[0].QDCsum[i];

            if( DEBUG ){
              printf("====================== event %u, event size %u\n", eventID, multi);
              printf("%6d, %12llu \n", event[0].energy, event[0].time);            
            }

            if( (rejectFlag & 0x8 ) && e[0] == 0 ) {
              event.clear();
              continue;
            }

            saveFile->cd();
            newtree->Fill();
            eventID ++;

            event.clear();
          }

          continue;

        }else{

          if( hitList[k].time - event.front().time <= timeWindow ){
            event.push_back(*(reader->data));
          }else{
            
            //save event
            if( DEBUG ) printf("====================== event %u, event size %lu\n", eventID, event.size());

            int nBGO = 0;
            int nClover = 0;
            int nGagg = 0;

            int count = 0 ;
            for( size_t p = 0; p < event.size(); p++ ) {
              if( (rejectFlag & 0x8 ) && event[p].energy == 0 ) continue;

              e[count]   = event[p].energy;
              e_t[count] = event[p].time;
              int index  = event[p].crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (event[p].slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + event[p].ch;
              id[count] = mapping[index];
              if( qdcFlag ) for( int i = 0; i < 8; i++) qdc[count][i] = event[p].QDCsum[i];

              if( DEBUG ) printf("%u | %3d, %6d, %12llu \n", count, id[count], e[count], e_t[count]);
              if(   0 <= id[count] && id[count] < 100 ) nClover ++;
              if( 100 <= id[count] && id[count] < 200 ) nBGO ++;
              if( 200 <= id[count] && id[count] < 400 ) nGagg ++;
              count ++;
            }

            multi = count;

            if( DEBUG ) printf("nBGO %d, nClover %d, nGagg %d \n", nBGO, nClover, nGagg);

            if( (rejectFlag & 0x1)  && nBGO > 0 ) {
              event.clear();
              event.push_back(*(reader->data));
              continue;
            }

            if( (rejectFlag & 0x2)  && nClover == 0 ) {
              event.clear();
              event.push_back(*(reader->data));
              continue;
            }

            if( (rejectFlag & 0x4)  && nGagg == 0 ) {
              event.clear();
              event.push_back(*(reader->data));
              continue;
            }

            if( multi > 0 ){
              if( DEBUG ) printf("---------------> fill tree\n");
              saveFile->cd();
              newtree->Fill();
              eventID ++;
            }

            //clear event
            event.clear();
            event.push_back(*(reader->data));

          }
        }
      
      }

    }while(true);

    tEnd = reader->data->time;
    delete reader;
  }

  //save the last event 
  if( timeWindow >= 0 ){
    multi = event.size();
    for( size_t p = 0; p < multi; p++ ) {
      if( DEBUG ) printf("%lu | %6d, %12llu \n", p, event[p].energy, event[p].time);
      int index  = event[p].crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (event[p].slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + event[p].ch;
      id[p] = mapping[index];
      e[p]   = event[p].energy;
      e_t[p] = event[p].time;
      if( qdcFlag ) for( int i = 0; i < 8; i++) qdc[p][i] = event[p].QDCsum[i];
    }
    saveFile->cd();
    newtree->Fill();
    eventID ++;
  }

  newtree->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 = %u by event builder (%llu in tree)\n", eventID , newtree->GetEntriesFast());
  double tDuration_sec = (tEnd - tStart) * 1e-9;
  printf("     first timestamp = %20llu ns\n", tStart);
  printf("      last timestamp = %20llu ns\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");

  saveFile->Close();


  return 0;
}


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;
}

unsigned long long getTime_ns(){
  std::chrono::high_resolution_clock::time_point currentTime = std::chrono::high_resolution_clock::now();
  std::chrono::nanoseconds nanoseconds = std::chrono::duration_cast<std::chrono::nanoseconds>(currentTime.time_since_epoch());
  return nanoseconds.count();
}

std::string trimSpaces(const std::string& str) {
  std::string trimmed = str;
  trimmed.erase(std::remove_if(trimmed.begin(), trimmed.end(), ::isspace), trimmed.end());
  return trimmed;
}

std::vector<std::string> split(const std::string& str, char delimiter) {
  std::vector<std::string> tokens;
  std::stringstream ss(str);
  std::string token;

  while (std::getline(ss, token, delimiter)) {
      tokens.push_back(trimSpaces(token));
  }

  return tokens;
}