XIAEventBuilder/xia2root.cpp

#include <stdio.h>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string.h>
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TMath.h"
#include <vector>

#define NUMDET 64  /// number of detector
#define STARTDETID 15

std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){

  std::vector<std::string> output;

  size_t pos;
  do{
    pos = tempLine.find(splitter); /// fine splitter
    if( pos == 0 ){ ///check if it is splitter again
      tempLine = tempLine.substr(pos+1);
      continue;
    }

    std::string secStr;
    if( pos == std::string::npos ){
      secStr = tempLine;
    }else{
      secStr = tempLine.substr(0, pos+shift);
      tempLine = tempLine.substr(pos+shift);
    }

    ///check if secStr is begin with space
    while( secStr.substr(0, 1) == " "){
      secStr = secStr.substr(1);
    };
    
    ///check if secStr is end with space
    while( secStr.back() == ' '){
      secStr = secStr.substr(0, secStr.size()-1);
    }

    output.push_back(secStr);
    //printf(" |%s---\n", secStr.c_str());
    
  }while(pos != std::string::npos );

  return output;
}

std::vector<std::vector<double>> LoadCorrectionParameters(TString corrFile){

  printf("==================== load correction parameters : %s", corrFile.Data());
  std::ifstream file;
  file.open(corrFile.Data());
  
  std::vector<std::vector<double>> corr;
  corr.clear();
  
  std::vector<double> detCorr;
  detCorr.clear();
  
  if( file.is_open() ){  
    
    while( file.good() ){
    
      std::string line;
      getline(file, line);
    
      if( line.substr(0,1) == "#" ) continue;
      if( line.substr(0,2) == "//" ) continue;
      if( line.size() == 0 ) continue;
      
      std::vector<std::string> temp = SplitStr(line, " ");
      
      detCorr.clear();
      for( int i = 0; i < (int) temp.size() ; i++){
        detCorr.push_back(std::stod(temp[i]));
      }
      corr.push_back(detCorr);
    }
    
    file.close();
    
    printf(".... done\n");
    printf("===== correction parameters \n");
    for( int i = 0; i < (int) corr.size(); i++){
      printf("det : %2d | ", i );
      int len = (int) corr[i].size();
      for( int j = 0; j < len - 1 ; j++){
        printf("%6.2f, ", corr[i][j]);
      }
      printf("%6.2f\n", corr[i][len-1]);
    }
  
  }else{
    printf(".... fail\n");
  }
  
  return corr;
}


//###################################################################################
//################                                        ###########################
//################        main                            ###########################
//################                                        ###########################
//###################################################################################
int main(int argn,char **argv) {
  if ( argn == 1 ) {
    printf("Usage: \n");
    printf("%s file_in.evt  raw_Opt  timeWidow   correctionFile\n", argv[0]);
    printf("                            |        |        |\n");
    printf("                            |        |        + correction file, row for det, col for order of correction\n");
    printf("                            |        |\n");
    printf("                            |        + when build event, event build window, 1 = 10 ns, default 100\n");
    printf("                            + default 0 = raw, 1 = event build \n");
    ///   std::cerr<<"Usage:\n "<<argv[0]<<" filein.evt fileout correctfilename STARTDETIDnumber maxtime start_buf# stop_buf#\n  
    ///    Converts physics buffers (type 30) from  \'file.evt\' to \'fileout.ev2\' from \'start_buf#\' until\n 
    ///   'stop_buf#'  or whole file if no event numbers are given or \'start buf#\' = 0\n  
    ///   STARTDETIDnum is the number of the dE ADC on a scale of 1 to n.\n  
    ///    maxtime defaults to 100 if not specified\n  
    ///   **Use \'NULL\' for the output file to avoid writing an ev2 file. \n
    ///   Note that the start and end buffer numbers are for the input where each adc \n
    ///   read forms another buffer.\n This version sets the time in event of the start detector to 100 and all other times relative to that. 
    ///   If a correct file is specified, it dorrects the detectors in the file for walk.  
    ///   It does not write out events which do not contain the start detector.\n
    ///   Also it only writes out events in which every ADC times in within the limits specified in the correct.tab file."<<std::endl;
    exit(0);
  }

  printf("=======================================================\n");
  printf("===             XIA evt file to CERN ROOT           ===\n");
  printf("=======================================================\n");
  printf("The start detector number is %d\n", STARTDETID);

  unsigned long long blockNum=0;  /// this is total block
  unsigned long long block30Num=0;  /// this is number of block-30
  unsigned long long eventID=0; /// this is eventID
  
  ///for( int i = 0 ; i < argn; i++) printf("............. %s \n", argv[i]);
  
  int rawOpt = 0;
  if ( argn >= 3 ) rawOpt = atoi(argv[2]);
  
  int timeWindow = 100;
  if ( argn >= 4 ) timeWindow = atoi(argv[3]);
  
  TString corrFileName = "";
  bool hasCorr = false;
  if ( argn == 5 ) {
    corrFileName = argv[4];
    hasCorr = true;
  }
  
  FILE *infile=fopen(argv[1],"r");
  if (infile==NULL) {
    printf("cannot open file : %s \n", argv[1]);
    ///std::cerr<<"Problem opening "<<argv[1]<<std::endl;
    exit(0);
  }

  TString inFileName = argv[1];
  TString outFileName = inFileName;
  outFileName.Remove(inFileName.First('.'));
  if( rawOpt == 0 ) outFileName.Append("_raw");
  outFileName.Append(".root");
  
  printf(" In file : %s \n",  inFileName.Data());
  printf("Out file : %s \n", outFileName.Data());

  std::vector<std::vector<double>> corr;
  if( hasCorr ) corr = LoadCorrectionParameters(corrFileName);

  int chan,slot,chann;
  int pu; /// pile up
  int energy;
  double cEnergy;
  unsigned long long evtime;
  unsigned short cfd;  

  int pileupcount = 0;
  int zerocount = 0;
  int PileUp[64];
  
  const unsigned long maskpu = 2147483648;
  const unsigned long multiplier = 4294967296LL;
  
  double energyA[NUMDET];
  double cEnergyA[NUMDET];
  unsigned long long  timeA[NUMDET];
  int puA[NUMDET];
  long long diffTimeA[NUMDET];
  unsigned short cfdA[NUMDET];
  int multi = 0; /// multipicilty in an event
  int detMulti[NUMDET]; /// multiplicity in a detector in an event

  TFile * outFile = new TFile(outFileName, "RECREATE");
  outFile->cd();
  TTree * tree = new TTree("tree", "tree");
    
  tree->Branch("eventID", &eventID, "event_number/l");

  if ( rawOpt == 0 ){ /// when save raw data
    tree->Branch("chan",     &chan, "chan/I");
    tree->Branch("slot",     &slot, "slot/I");
    tree->Branch("chann",   &chann, "channel number/I");
    tree->Branch("pu",         &pu, "pile-up/I");
    tree->Branch("energy", &energy, "energy/I");
    if( hasCorr) tree->Branch("cEnergy", &cEnergy, "corrected_energy/D");
    tree->Branch("time",   &evtime, "timestamp/l");
    tree->Branch("cfd",       &cfd, "cfd/s");
  }else{ /// when build event by time-window
    tree->Branch("energy",    energyA, Form("energy[%d]/D", NUMDET));
    if( hasCorr) tree->Branch("cEnergy",    cEnergyA, Form("corrected_energy[%d]/D", NUMDET));
    tree->Branch("time",        timeA, Form("timestamp[%d]/l", NUMDET));
    tree->Branch("dtime",   diffTimeA, Form("diff_time[%d]/L", NUMDET));
    tree->Branch("pu",            puA, Form("pile_up[%d]/I", NUMDET));
    tree->Branch("cfd",          cfdA, Form("cfd[%d]/I", NUMDET));
    tree->Branch("multi",      &multi, "multiplicity/I");
    tree->Branch("detMulti", detMulti, Form("det_multiplicity[%d]/I", NUMDET));
  }

  ///change this for 64bit compiler long *bufsam=NULL;
  
  //clear energy and time array
  for( int i = 0; i < NUMDET; i++){
    energyA[i] = TMath::QuietNaN();
    cEnergyA[i] = TMath::QuietNaN();
    timeA[i] = 0;
    diffTimeA[i] = -999;
    cfdA[i] = 0;
    puA[i] = -1;
    detMulti[i] = 0;
  }
  multi = 0;
  
  unsigned long long startTime = 0;
  long long diffTime = 0;

  int bread = 1; 
  int bsam = 2;
  long * bufsiz=new long[bsam];
  unsigned int *bufsam = NULL;

  printf("============ Start looping events | build event ? %s", rawOpt == 0 ? "No" : "Yes");
  if( rawOpt == 1 ) {
    printf(" time window : +- %d click\n", timeWindow);
  }else{
    printf("\n");
  }
  while ( !feof(infile) ) {
  
    // get buffer size
    ///change long -> int for 64 bit
    fread(bufsiz,sizeof(int),bread,infile); /// read 1 int (4 byte) from infile and save to bufsize
    int bsize = bufsiz[0] -4 ;
    if (feof(infile)) break;
    blockNum ++;

    ///change for 64bit  bufsam=new long[bsize/4+1];
    bufsam = new unsigned int[bsize/4+1];
    fread((char*)bufsam, 1, bsize, infile); /// read bsize of 1 byte from infile and save to char
    ///printf("============ bsize : %d \n", bsize);
    ///for( int i = 0; i < bsize; i++) printf("%d, ", bufsam[i]);
    ///printf("\n");

    if (bufsam[0] == 30) {
      block30Num ++;
      chan = (bufsam[2]) & (15);
      slot = ((bufsam[2]) & (240))/16;
      chann = (slot - 2)*16 + chan + 1;
      pu = ((bufsam[2]) & (maskpu))/maskpu;
      energy = ((bufsam[5]) & 65535);
      unsigned long long evtimehi = ((bufsam[4]) & 65535);
      unsigned long long evtimelo = bufsam[3];
      evtime = evtimelo + multiplier*evtimehi;
      cfd = bufsam[4]/65536;

      if ( energy == 0 ) zerocount++;
      if ( pu > 0 ) pileupcount++;
      if ((pu > 0 ) && ( chann > 0 ) && ( chann < 65 )) PileUp[chann-1]++;

      if( blockNum % 100000 == 0 ) printf(".");
      ///if( blockNum % 100000 == 0 ) printf("%llu \n", blockNum);
      ///if( block30Num < 50) printf("b30: %10llu, chan: %d, slot: %d, chann: %2d, pu: %2d, energy: %5d, evtime: %13llu, cfd: %d\n", block30Num, chan, slot, chann, pu, energy, evtime, cfd);

      /// energy correction
      if ( hasCorr ){
        cEnergy = 0;
        int order = (int) corr[chann-1].size();
        for( int i = 0; i < order ; i++){
          cEnergy += corr[chann-1][i] * TMath::Power((double)energy, i);
        }
      } 

      if ( rawOpt == 0 ) {
        eventID++;
        outFile->cd();
        tree->Fill();
      }else{ /// build event

        if ( startTime == 0 ) startTime = evtime;
      
        diffTime = evtime - startTime;
        
        if( -timeWindow < diffTime && diffTime < timeWindow ){

          if( !TMath::IsNaN(energyA[chann-1]) ) detMulti[chann-1] ++;
          energyA[chann-1] = energy;
          cEnergyA[chann-1] = cEnergy;
          timeA[chann-1] = evtime;
          diffTimeA[chann-1] = diffTime;
          puA[chann-1] = pu;
          detMulti[chann-1]++;
          multi++;
          
          
        }else{
          /// fill tree 
          eventID++;
          outFile->cd();
          tree->Fill();
          ///  clear energy and time array
          multi = 0;
          for( int i = 0; i < NUMDET; i++){
            energyA[i] = TMath::QuietNaN();
            cEnergyA[i] = TMath::QuietNaN();
            timeA[i] = 0;
            diffTimeA[i] = -999;
            puA[i] = -1;
            detMulti[i] = 0;
            cfdA[i] = 0;
          }
          
          /// fill the 1st data of a new event
          startTime = evtime;
          energyA[chann-1] = energy;
          cEnergyA[chann-1] = cEnergy;
          timeA[chann-1] = evtime;
          diffTimeA[chann-1] = 0;
          puA[chann-1] = pu;
          detMulti[chann-1]++;
          multi++;
        }
      }
      
    } ///end if bufsam[0]=30
    
  } 

  delete [] bufsiz;
  delete [] bufsam;
  fclose(infile);

  printf("\n============ end of event loop, totoal block read: %llu \n", blockNum);

  eventID++;
  outFile->cd();
  tree->Write();
  outFile->Close();
  
  
  //========================= Print summary
  printf("============================================\n");
  ///printf("         number of block: %llu\n", blockNum); 
  printf(" number of type 30 block: %llu\n", block30Num); 
  printf("             event built: %llu\n", eventID); 
  printf("============================================\n");
    
  return 0;
}