362 lines
12 KiB
C++
362 lines
12 KiB
C++
/*==================
|
|
|
|
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] [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(" 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());
|
|
}
|
|
|
|
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};
|
|
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("File-%ld %10lu block %10lu / %lu | %u \n", i, 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;
|
|
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( 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 ) {
|
|
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( 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;
|
|
}
|
|
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;
|
|
} |