diff --git a/armory/pxi-fsu-time-order.cpp b/armory/pxi-fsu-time-order.cpp index dacdfdd..f8ea24a 100644 --- a/armory/pxi-fsu-time-order.cpp +++ b/armory/pxi-fsu-time-order.cpp @@ -44,22 +44,67 @@ #define M2_SUB_EVENTS 5000 #define M3_SUB_EVENTS 20000 #define M4_SUB_EVENTS 50000 + #define M5_SUB_EVENTS 100000 #define MAX_SUB_EVENTS 200000 #define MAX_MALLOC 4000*1024*1024L //2GB +#define MAXLONGLONGINT 9223372036854775807 + +//TODO change to load file +#include "../mapping.h" + struct subevent { long long int timestamp; int length; //unit = words with 4 bytes per word + int detID; unsigned int *data; }; struct subevent *subevents[MAX_ID]; int nevts[MAX_ID], iptr[MAX_ID]; int maxevts[MAX_ID]; - + + +void swap64(long long int * a, long long int *b){ + long long int t = *a; + *a = *b; + *b = t; +} + +void swapInt(int * a, int *b){ + int t = *a; + *a = *b; + *b = t; +} + +int partition(long long int arr[], int index [], int low, int high){ + + long long int pivot = arr[high]; + int i = (low -1); + + for(int j = low; j <= high -1 ; j++){ + if( arr[j] < pivot ){ + i++; + swap64(&arr[i], &arr[j]); + swapInt(&index[i], &index[j]); + } + } + swap64(&arr[i+1], &arr[high]); + swapInt(&index[i+1], &index[high]); + return i+1; +} + +void quickSort(long long int arr[], int index[], int low, int high){ + if (low < high){ + int pi = partition(arr, index, low, high); + quickSort(arr, index, low, pi - 1); + quickSort(arr, index, pi + 1, high); + } +} + int main(int argc, char **argv) { @@ -73,6 +118,7 @@ int main(int argc, char **argv) { int pause=0; long long int nwords=0, evts_tot_read=0, evts_tot_write=0; + long long int evts_tot_drop = 0; long long int time=0,time_old=0; int length=0; @@ -91,11 +137,13 @@ int main(int argc, char **argv) { int min_id = -1; memset(nevts, 0, sizeof(nevts)); - memset(iptr, 0, sizeof(iptr)); + memset(iptr, 0, sizeof(iptr)); /// index of time + + long long int timeIndex[MAX_ID]; + int index[MAX_ID]; + int fillSize = 0; int i=0, j=0; - div_t e_div; - //open input event file if ((fpr = fopen(argv[1], "r")) == NULL) { @@ -114,6 +162,14 @@ int main(int argc, char **argv) { return 1; } + printf("open : \033[1;31m%s\033[m\n", argv[1]); + printf("save : \033[1;34m%s\033[m\n", filenameto); + + int eventWindow = 0; + + if( argc >= 3 ) eventWindow = atoi(argv[2]); + printf(" event build time window : %d ticks = %d ns \n", eventWindow, 10* eventWindow); + //check for lockfile, active PID, and event file for auto "online" mode detection FILE *FPLOCK; @@ -142,19 +198,15 @@ int main(int argc, char **argv) { fgets(line, 1024, FPPATH); //skip first line fgets(line, 1024, FPPATH); //need second line sscanf(line,"%s\n", onlinefile); - fclose(FPPATH); - if (filename == NULL) { - if (strcmp(onlinefile,argv[1]) == 0) { - online = 1; - } - }else { - if (strcmp(onlinefile,filename+1) == 0) { - online = 1; - } - } + fclose(FPPATH); + if (filename == NULL) { + if (strcmp(onlinefile,argv[1]) == 0) online = 1; + }else { + if (strcmp(onlinefile,filename+1) == 0) online = 1; + } } - } - } + } + } if (online == 1) printf("Auto Mode: \x1B[32mOnline\x1B[0m\n"); else printf("Auto Mode: \x1B[32mOffline\x1B[0m\n"); @@ -164,11 +216,8 @@ int main(int argc, char **argv) { fprsize = fprsize_orig = ftell(fpr); fseek(fpr, fprpos, SEEK_SET); - - - //Get memory for default number of subevents per channel id - for (i=0; i= 0) { - FPPATH = fopen(pathfile, "r"); - if (FPPATH != NULL) { - fgets(line, 1024, FPPATH); //skip first line - fgets(line, 1024, FPPATH); //need second line - sscanf(line,"%s\n", onlinefile); - fclose(FPPATH); - if (filename == NULL) { - if (strcmp(onlinefile,argv[1]) == 0) online = 1; - } - else { - if (strcmp(onlinefile,filename+1) == 0) online = 1; - } - } - } - } - } //end auto online + //check for lock file and active PID + FPLOCK = fopen(lockfile, "r"); + if (FPLOCK != NULL) { + fscanf(FPLOCK, "%d", &lockpid); + fclose(FPLOCK); + if (getpgid(lockpid) >= 0) { + FPPATH = fopen(pathfile, "r"); + if (FPPATH != NULL) { + fgets(line, 1024, FPPATH); //skip first line + fgets(line, 1024, FPPATH); //need second line + sscanf(line,"%s\n", onlinefile); + fclose(FPPATH); + if (filename == NULL) { + if (strcmp(onlinefile,argv[1]) == 0) online = 1; + }else { + if (strcmp(onlinefile,filename+1) == 0) online = 1; + } + } + } + } + } //end auto online - //read 4-byte header - if (fread(subhead, sizeof(subhead), 1, fpr) != 1) break; - nwords = nwords + HEADER_LENGTH; - chn = subhead[0] & 0xF; - sln = (subhead[0] & 0xF0) >> 4; - crn = (subhead[0] & 0xF00) >> 8; - id = crn*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (sln-BOARD_START)*MAX_CHANNELS_PER_BOARD + chn; - length = (subhead[0] & 0x7FFE0000) >> 17; //unit = words with 4 bytes per word - time = ( (long long int)(subhead[2] & 0xFFFF) << 32) + subhead[1]; - if (id > idmax) idmax = id; + //read 4-byte header + if (fread(subhead, sizeof(subhead), 1, fpr) != 1) break; + nwords = nwords + HEADER_LENGTH; + chn = subhead[0] & 0xF; + sln = (subhead[0] & 0xF0) >> 4; + crn = (subhead[0] & 0xF00) >> 8; + id = crn*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (sln-BOARD_START)*MAX_CHANNELS_PER_BOARD + chn; + length = (subhead[0] & 0x7FFE0000) >> 17; //unit = words with 4 bytes per word + time = ( (long long int)(subhead[2] & 0xFFFF) << 32) + subhead[1]; + if (id > idmax) idmax = id; } @@ -251,54 +302,57 @@ int main(int argc, char **argv) { //Expand memory for more events (careful when final is to left of initial in circular buffer) if ( maxevts[id] - nevts[id] == 1 && totmem < MAX_MALLOC) { - - if (maxevts[id] == DEF_SUB_EVENTS) {evts_old = DEF_SUB_EVENTS; evts_new = M1_SUB_EVENTS;} - if (maxevts[id] == M1_SUB_EVENTS) {evts_old = M1_SUB_EVENTS; evts_new = M2_SUB_EVENTS;} - if (maxevts[id] == M2_SUB_EVENTS) {evts_old = M2_SUB_EVENTS; evts_new = M3_SUB_EVENTS;} - if (maxevts[id] == M3_SUB_EVENTS) {evts_old = M3_SUB_EVENTS; evts_new = M4_SUB_EVENTS;} - if (maxevts[id] == M4_SUB_EVENTS) {evts_old = M4_SUB_EVENTS; evts_new = M5_SUB_EVENTS;} - if (maxevts[id] == M5_SUB_EVENTS) {evts_old = M5_SUB_EVENTS; evts_new = MAX_SUB_EVENTS;} + + if (maxevts[id] == DEF_SUB_EVENTS) {evts_old = DEF_SUB_EVENTS; evts_new = M1_SUB_EVENTS;} + if (maxevts[id] == M1_SUB_EVENTS) {evts_old = M1_SUB_EVENTS; evts_new = M2_SUB_EVENTS;} + if (maxevts[id] == M2_SUB_EVENTS) {evts_old = M2_SUB_EVENTS; evts_new = M3_SUB_EVENTS;} + if (maxevts[id] == M3_SUB_EVENTS) {evts_old = M3_SUB_EVENTS; evts_new = M4_SUB_EVENTS;} + if (maxevts[id] == M4_SUB_EVENTS) {evts_old = M4_SUB_EVENTS; evts_new = M5_SUB_EVENTS;} + if (maxevts[id] == M5_SUB_EVENTS) {evts_old = M5_SUB_EVENTS; evts_new = MAX_SUB_EVENTS;} - if (maxevts[id]==evts_old && totmem + (evts_new-evts_old)*(sizeof(struct subevent) + sizeof(unsigned int)*length) < MAX_MALLOC) { - subevents[id] = (struct subevent *) realloc(subevents[id], sizeof(struct subevent)*evts_new); - if (subevents[id] == NULL) { - printf("realloc failed\n"); - return -1; - } - totmem = totmem - sizeof(struct subevent)*evts_old + sizeof(struct subevent)*evts_new; - maxevts[id] = evts_new; - for (j=evts_old; j evts_old) { - for (j=0; j evts_old) { + for (j=0; j= maxevts[id]) j -= maxevts[id]; + + subevents[id][j].timestamp = time; + subevents[id][j].detID = mapping[id]; + if (subevents[id][j].data == NULL) { + subevents[id][j].data = (unsigned int *) malloc(sizeof(unsigned int)*length); + if (subevents[id][j].data == NULL) { + printf("malloc failed\n"); + return -1; + } + totmem += sizeof(unsigned int)*length; + }else if (length != subevents[id][j].length) { //not needed anymore since always free data after use now. Keep for future ... + subevents[id][j].data = (unsigned int *) realloc(subevents[id][j].data, sizeof(unsigned int)*length); + if (subevents[id][j].data == NULL) { + printf("realloc failed\n"); + return -1; + } + totmem = totmem - sizeof(unsigned int)*subevents[id][j].length + sizeof(unsigned int)*length; + } - j = nevts[id] + iptr[id]; - if (j >= maxevts[id]) j -= maxevts[id]; - - subevents[id][j].timestamp = time; + subevents[id][j].length = length; - if (subevents[id][j].data == NULL) { - subevents[id][j].data = (unsigned int *) malloc(sizeof(unsigned int)*length); - if (subevents[id][j].data == NULL) { - printf("malloc failed\n"); - return -1; - } - totmem += sizeof(unsigned int)*length; - } - else if (length != subevents[id][j].length) { //not needed anymore since always free data after use now. Keep for future ... - subevents[id][j].data = (unsigned int *) realloc(subevents[id][j].data, sizeof(unsigned int)*length); - if (subevents[id][j].data == NULL) { - printf("realloc failed\n"); - return -1; - } - totmem = totmem - sizeof(unsigned int)*subevents[id][j].length + sizeof(unsigned int)*length; - } + if (length>HEADER_LENGTH) { + if (fread(subevents[id][j].data + HEADER_LENGTH, (length-HEADER_LENGTH)*sizeof(int), 1, fpr) != 1) break; + nwords = nwords + (length-HEADER_LENGTH); + } + + for (i=0; i < HEADER_LENGTH; i++) { + subevents[id][j].data[i] = subhead[i]; + } - subevents[id][j].length = length; + nevts[id]++; + evts_tot_read++; - if (length>HEADER_LENGTH) { - if (fread(subevents[id][j].data + HEADER_LENGTH, (length-HEADER_LENGTH)*sizeof(int), 1, fpr) != 1) break; - nwords = nwords + (length-HEADER_LENGTH); - } - - for (i=0; i < HEADER_LENGTH; i++) { - subevents[id][j].data[i] = subhead[i]; - } - - nevts[id]++; - evts_tot_read++; - - } - else { - pause = 1; - break; + }else { + pause = 1; + break; } - } // end while for fill buffers + } // end while for fill buffers maxevts[id] ///////// - - ///////// - // write event with minimum time to file - timemin_old = timemin; - timemin = -1; - for (i=0; i < idmax + 1; i++) { //could be MAX_ID but limit ourselves to current max, idmax - if (nevts[i] > 0) { - if (timemin == -1) { - timemin = subevents[i][iptr[i]].timestamp; - time_old = timemin; - min_id = i; - } - else if (subevents[i][iptr[i]].timestamp < timemin) { - timemin = subevents[i][iptr[i]].timestamp; - time_old = timemin; - min_id = i; - } - } - } - - - if (timemin > -1) { - if (timemin < timemin_old) { - printf("\nWarning!!! timemin = %lld and timemin_old = %lld and min_id = %d\n", timemin, timemin_old, min_id); - outoforder++; - } - if (subevents[min_id][iptr[min_id]].data == NULL) {printf("Error: data = NULL\n"); return -1;} - fwrite(subevents[min_id][iptr[min_id]].data, sizeof(unsigned int)*subevents[min_id][iptr[min_id]].length, 1, fpw); - - //free data memory up until it's needed again - free(subevents[min_id][iptr[min_id]].data); - subevents[min_id][iptr[min_id]].data = NULL; - totmem -= sizeof(unsigned int)*subevents[min_id][iptr[min_id]].length; - subevents[min_id][iptr[min_id]].length = 0; - subevents[min_id][iptr[min_id]].timestamp = 0; + + //######################## FSU + // find group of event within timewindow, check is contain gamma. if not, throw away all. + + if( eventWindow > 0 ){ + + //quick sort of subevents[i][iptr[i]].timestamp, i = 0 , idmax +1 + for( i = 0 ; i < idmax + 1; i++) { - nevts[min_id]--; - if (++iptr[min_id] >= maxevts[min_id]) iptr[min_id] -= maxevts[min_id]; - evts_tot_write++; + if( nevts[i] > 0 ){ + timeIndex[i] = subevents[i][iptr[i]].timestamp; + }else{ + timeIndex[i] = MAXLONGLONGINT; + } + index[i] = i; + //printf("%3d , %llu, %d \n", i, timeIndex[i], index[i] ); + } + quickSort(timeIndex, index, 0, idmax); + for( i = 0 ; i < idmax + 1; i++) { + //printf("%3d , %llu , %d\n", i, timeIndex[i], index[i]); + } + + //reduce the index size + fillSize = 1; + for( i = 1; i < idmax + 1; i++){ + if( timeIndex[i] - timeIndex[0] < eventWindow) fillSize ++; + } + + + // display + if ( count < debugCount) { + //if ( count < debugCount || timeIndex[10] == MAXLONGLONGINT) { + printf("===============================================\n"); + + for( int i = 0; i < idmax+1; i++) { + printf("%3d, %llu, %d \n", i, timeIndex[i], index[i]); + if( i == fillSize - 1 ) printf("------------------- %d \n", fillSize); + } + + } + + if( timeIndex[0] == MAXLONGLONGINT ) break; + + //CHeck if fill evt.to for the data. + bool fillFlag = false; + + for( i = 0 ; i < fillSize; i++ ){ + if( count < debugCount ) printf("********************* %llu , detID : %d\n", timeIndex[0], subevents[index[i]][iptr[index[i]]].detID); + if( subevents[index[i]][iptr[index[i]]].detID < 100 ) fillFlag = true; + } + if( count < debugCount ) printf("=============== fillFlag : %d \n", fillFlag); + + for( i = 0; i < fillSize; i++){ + min_id = index[i]; + + if( fillFlag ){ + fwrite(subevents[min_id][iptr[min_id]].data, sizeof(unsigned int)*subevents[min_id][iptr[min_id]].length, 1, fpw); + + if( count < debugCount ) printf("filling subevents[%d][iprt[%d]] \n", min_id, min_id); + evts_tot_write++; + }else{ + evts_tot_drop++; + } + + + //free data memory up until it's needed again + if( count < debugCount ) printf(" Free subevents[%d][iprt[%d]] \n", min_id, min_id); + + free(subevents[min_id][iptr[min_id]].data); + subevents[min_id][iptr[min_id]].data = NULL; + totmem -= sizeof(unsigned int)*subevents[min_id][iptr[min_id]].length; + subevents[min_id][iptr[min_id]].length = 0; + subevents[min_id][iptr[min_id]].timestamp = 0; + subevents[min_id][iptr[min_id]].detID = -1; + + nevts[min_id]--; + if (++iptr[min_id] >= maxevts[min_id]) iptr[min_id] -= maxevts[min_id]; + + } + + count ++; + + }else{ + + ///////// + // write event with minimum time to file + timemin_old = timemin; + timemin = -1; + for (i=0; i < idmax + 1; i++) { //could be MAX_ID but limit ourselves to current max, idmax + if (nevts[i] > 0) { + if (timemin == -1) { + timemin = subevents[i][iptr[i]].timestamp; + time_old = timemin; + min_id = i; + }else if (subevents[i][iptr[i]].timestamp < timemin) { + timemin = subevents[i][iptr[i]].timestamp; + time_old = timemin; + min_id = i; + } + } + } + + if (timemin > -1) { + if (timemin < timemin_old) { + printf("\nWarning!!! timemin = %lld and timemin_old = %lld and min_id = %d\n", timemin, timemin_old, min_id); + outoforder++; + } + if (subevents[min_id][iptr[min_id]].data == NULL) {printf("Error: data = NULL\n"); return -1;} + + fwrite(subevents[min_id][iptr[min_id]].data, sizeof(unsigned int)*subevents[min_id][iptr[min_id]].length, 1, fpw); + + //free data memory up until it's needed again + free(subevents[min_id][iptr[min_id]].data); + subevents[min_id][iptr[min_id]].data = NULL; + totmem -= sizeof(unsigned int)*subevents[min_id][iptr[min_id]].length; + subevents[min_id][iptr[min_id]].length = 0; + subevents[min_id][iptr[min_id]].timestamp = 0; + + nevts[min_id]--; + if (++iptr[min_id] >= maxevts[min_id]) iptr[min_id] -= maxevts[min_id]; + evts_tot_write++; + + }else break; + ///////// + } - else break; - ///////// - //print statistics - //e_div=div(evts_tot_read,10000); - //if ( e_div.rem == 0) if( evts_tot_read % 10000 == 0 ) - printf("Malloc (%d MB) : evts in (\x1B[34m%lld\x1B[0m) : evts out (\x1B[32m%lld\x1B[0m) : diff (\x1B[31m%lld\x1B[0m)\r", (totmem)/1024/1024, evts_tot_read, evts_tot_write, evts_tot_read-evts_tot_write); + printf("Malloc (%d MB) : evts in (\x1B[34m%lld\x1B[0m) : evts out (\x1B[32m%lld\x1B[0m) : evts drop (\x1B[32m%lld\x1B[0m) : diff (\x1B[31m%lld\x1B[0m)\r", + (totmem)/1024/1024, evts_tot_read, evts_tot_write, evts_tot_drop, evts_tot_read-evts_tot_write - evts_tot_drop); } //end main while diff --git a/armory/to2root.cpp b/armory/to2root.cpp index 38bdd79..5d67a7e 100644 --- a/armory/to2root.cpp +++ b/armory/to2root.cpp @@ -102,13 +102,14 @@ int main(int argc, char **argv) { int countGP = 0; //gamma-particle coincident double totalDataSize = 0; + int outFileCount = 0; for( int i = 0; i < nFile; i++){ evt->OpenFile(inFileName[i]); if( evt->IsOpen() == false ) continue; - printf("==================================================== %d / %d\n", i, nFile); + printf("==================================================== %d / %d\n", i+1, nFile); printf("\033[1;31m%s \033[m\n", inFileName[i].Data()); int pos = inFileName[i].Last('/'); @@ -141,13 +142,13 @@ int main(int argc, char **argv) { if (tdif > timeWindow) { //Gate - if( multiCry > 0 && multiGagg > 0 ) { - - outRootFile->cd(); - tree->Fill(); - - countGP++; - } + //if( multiCry > 2 && multiGagg ==0 ) { + // + outRootFile->cd(); + tree->Fill(); + // + // countGP++; + //} evID ++; @@ -200,8 +201,38 @@ int main(int argc, char **argv) { printf(" ----------- root file size : %.3f GB\n", rootFileSize); printf(" ---------- total read size : %.3f GB\n", totalDataSize); printf(" ----------- reduction rate : %.3f %%\n", rootFileSize*100./totalDataSize); - if( rootFileSize > 3.0 ) break; - + + if( rootFileSize > 3.0 ) { + break; + } + + ///try to open a new root file when file size > 2 GB + /*if( rootFileSize > 2.0 ) { + + outRootFile->Close(); + delete outRootFile; + delete tree; + outFileCount += 1; + + if( outFileCount > 5 ) break; + + TString outFileName2 = outFileName; + outFileName2.Insert(outFileName.Sizeof() - 6, Form("_%03d",outFileCount)); + outRootFile = new TFile( outFileName2, "recreate"); + + tree = new TTree("tree", "tree"); + tree->Branch("evID", &evID, "event_ID/l"); + tree->Branch("multi", &multi, "multi/I"); + tree->Branch("detID", id, "detID[multi]/I"); + tree->Branch("e", e, "e[multi]/D"); + tree->Branch("e_t", e_t, "e_timestamp[multi]/l"); + tree->Branch("qdc", qdc, "qdc[multi][8]/I"); + tree->Branch("multiCry", &multiCry, "multiplicity_crystal/I"); + tree->Branch("multiGagg", &multiGagg, "multiplicity_GAGG/I"); + tree->Branch("runID", &runID, "runID/I"); + + }*/ + }