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 "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TBenchmark.h"
#include "TStopwatch.h"
#include "TTreeIndex.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 <chrono>
// inline 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();
// }

#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] [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, 1 = reject no gamma, 2 = reject no GAGG. see mapping.h\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]);
  TString outFileName = argv[4];

  std::vector<std::string> inFileList;
  for( int i = 5; 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());
  }

  printf("================================== Digesting Mapping file (to be impletmented\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};
  newtree->Branch("id",       id, "id[multi]/i" );
  newtree->Branch("e",         e, "e[multi]/I" );
  newtree->Branch("e_t",     e_t, "e_timestamp[multi]/l");

  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();
      printf(" %lu block %lu / %lu | %u \n", hitList.size(), blockCount, totalBlock, eventID);
      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;
            e[0]   = event[0].energy;
            e_t[0] = event[0].time;

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

            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
            multi = event.size();

            if( DEBUG ) printf("====================== event %u, event size %u\n", eventID, multi);
            for( size_t p = 0; p < multi; p++ ) {
              if( DEBUG ) printf("%lu | %6d, %12llu \n", p, event[p].energy, event[p].time);
              e[p]   = event[p].energy;
              e_t[p] = event[p].time;
            }

            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);
      e[p]   = event[p].energy;
      e_t[p] = event[p].time;
    }
    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;
}