/*================== Thie event builder both sort and build event. skip the *.to file. ===================*/ #include #include #include #include #include #include #include #include /** struct timeval, select() */ #include #include #include // for std::remove_if #include // 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 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 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 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 hitList; evtReader * reader = nullptr; std::vector 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(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 split(const std::string& str, char delimiter) { std::vector tokens; std::stringstream ss(str); std::string token; while (std::getline(ss, token, delimiter)) { tokens.push_back(trimSpaces(token)); } return tokens; }