diff --git a/.gitignore b/.gitignore index 0475b57..9a20272 100644 --- a/.gitignore +++ b/.gitignore @@ -25,7 +25,6 @@ ti74pt7a fsu_run_2021 data Data -obsolete raw_data root_data diff --git a/armory/obsolete/EventBuilder.cpp b/armory/obsolete/EventBuilder.cpp new file mode 100644 index 0000000..c174946 --- /dev/null +++ b/armory/obsolete/EventBuilder.cpp @@ -0,0 +1,180 @@ +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TTree.h" +#include "TMath.h" +#include "TBenchmark.h" +#include "TStopwatch.h" +#include "TTreeIndex.h" + +#include "../mapping.h" + +#define MAXMULTI 100 + +int main(int argn, char **argv){ + printf("=====================================\n"); + printf("=== Event Builder from *_raw.root ===\n"); + printf("=====================================\n"); + + if (argn != 2 && argn != 3 && argn != 4 ) { + printf("Usage :\n"); + printf("%s [_raw.root File] \n", argv[0]); + printf(" timeWindows : default = 100 \n"); + printf(" SaveFileName : default is *.root \n"); + return 1; + } + + TString inFileName = argv[1]; // need to check name + int timeWindow = 100; + if( argn >= 3 ) timeWindow = atoi(argv[2]); + + printf(">>> Opening input %s \n", inFileName.Data()); + TFile * inFile = new TFile(inFileName, "READ"); + if( inFile->IsOpen() == false ) { + printf("!!!! cannot open file %s \n", inFileName.Data()); + return 0; + } + + TTree * tree = (TTree *) inFile->Get("tree"); + + Long64_t evID; + UShort_t detID; + UShort_t energy; + ULong64_t energy_t; + + TBranch *b_data_ID; //! + TBranch *b_ID; //! + TBranch *b_energy; //! + TBranch *b_energy_timestamp; //! + + tree->SetBranchAddress("evID", &evID, &b_data_ID); + tree->SetBranchAddress("id", &detID, &b_ID); + tree->SetBranchAddress("e", &energy, &b_energy); + tree->SetBranchAddress("e_t", &energy_t, &b_energy_timestamp); + + Long64_t totnumEntry = tree->GetEntries(); + + printf(" total Entry : %lld \n", totnumEntry); + + printf(" event Build window: %d tick = %d nsec \n", timeWindow, timeWindow * 10); + + printf(">>> Buidling Index using the timestamp\n"); + tree->BuildIndex("e_t"); + TTreeIndex *in = (TTreeIndex*) tree->GetTreeIndex(); + Long64_t * index = in->GetIndex(); + + ULong64_t time0; //time-0 for each event + int timeDiff; + + TString outFileName = inFileName; + outFileName.Remove(inFileName.First("_raw")); + outFileName.Append(".root"); + if( argn >=4 ) outFileName = argv[3]; + + printf(">>> out File name : \033[1;31m%s\033[m\n", outFileName.Data()); + + printf(">>> Create output tree\n"); + TFile * saveFile = new TFile(outFileName, "recreate"); + saveFile->cd(); + TTree * newtree = new TTree("tree", outFileName); + + Int_t eventID = 0 ; + Int_t multi = 0; /// this is total multipicilty for all detectors + newtree->Branch("multi", &multi, "multi/I"); + + + newtree->Branch("evID", &eventID, "event_ID/l"); + + Int_t multiCry = 0 ; /// thi is total multiplicity for all crystal + newtree->Branch("multiCry", &multiCry, "multiplicity_crystal/I"); + + int id[MAXMULTI] = {0}; + double e[MAXMULTI] = {TMath::QuietNaN()}; + ULong64_t e_t[MAXMULTI] = {0}; + newtree->Branch("id", id, "id[multi]/I" ); + newtree->Branch("e", e, "e[multi]/D" ); + newtree->Branch("e_t", e_t, "e_timestamp[multi]/l"); + + printf("================== Start processing....\n"); + Float_t Frac = 0.05; ///Progress bar + TStopwatch StpWatch; + StpWatch.Start(); + + int multiOverflow = 0; + + for( Long64_t entry = 0; entry < totnumEntry; entry++){ + + /*********** Progress Bar ******************************************/ + if (entry>totnumEntry*Frac-1) { + TString msg; msg.Form("%llu", totnumEntry/1000); + int len = msg.Sizeof(); + printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n", + Frac*100, len, entry/1000,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); + StpWatch.Start(kFALSE); + Frac+=0.1; + } + + Long64_t ev = index[entry]; + + b_ID->GetEntry(ev, 0); + b_energy->GetEntry(ev, 0); + b_energy_timestamp->GetEntry(ev, 0); + + if( time0 == 0) { + time0 = energy_t; + multi = 0; + } + timeDiff = (int) (energy_t - time0); + + if( timeDiff < timeWindow ) { + + if( multi > MAXMULTI ){ + multiOverflow++; + }else{ + id[multi] = detID; + e[multi] = energy; + e_t[multi] = energy_t; + multi ++; + if( detID < NCRYSTAL ) multiCry++; + } + }else{ + ///---- end of event + saveFile->cd(); + newtree->Fill(); + eventID ++; + + ///---- clear data + for( int i = 0; i < MAXMULTI ; i ++){ + id[i] = 0; + e[i] = TMath::QuietNaN(); + e_t[i] = 0; + } + multi = 0; + multiCry = 0; + + + /// fill 1st data of an event + time0 = energy_t; + timeDiff = 0; + + id[multi] = detID; + e[multi] = energy; + e_t[multi] = energy_t; + multi++; + if( detID < NCRYSTAL ) multiCry++; + } + } + + printf("============================== finished.\n"); + + saveFile->cd(); + newtree->Write(); + saveFile->Close(); + + printf(" total number of event Built : %d \n", eventID); + printf(" total event has multi > %6d : %d \n", MAXMULTI, multiOverflow); +} diff --git a/armory/obsolete/EventBuilder_seperated.cpp b/armory/obsolete/EventBuilder_seperated.cpp new file mode 100644 index 0000000..45fd3cf --- /dev/null +++ b/armory/obsolete/EventBuilder_seperated.cpp @@ -0,0 +1,206 @@ +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TTree.h" +#include "TMath.h" +#include "TBenchmark.h" +#include "TStopwatch.h" +#include "TTreeIndex.h" + +#include "../mapping.h" + +Int_t eventID = 0 ; +double e[NCRYSTAL]; +ULong64_t e_t[NCRYSTAL]; +double bgo[NBGO]; +ULong64_t bgo_t[NBGO]; +Short_t other[NOTHER]; +Short_t multi; + +void ClearTreeData(){ + + for( int i = 0; i < NCRYSTAL; i++){ + e[i] = TMath::QuietNaN(); + e_t[i] = 0; + //pileup[i] = 0; + //hit[i] = 0; + } + for( int i = 0; i < NBGO; i++) { + bgo[i] = TMath::QuietNaN(); + bgo_t[i] = 0 ; + } + for( int i = 0; i < NOTHER; i++) { + other[i] = TMath::QuietNaN(); + } + multi = 0; +} + +int main(int argn, char **argv){ + printf("=====================================\n"); + printf("=== Event Builder ===\n"); + printf("=====================================\n"); + + if (argn != 2 && argn != 3 && argn != 4 ) { + printf("Usage :\n"); + printf("%s [_raw.root File] \n", argv[0]); + printf(" timeWindows : default = 100 \n"); + printf(" SaveFileName : default is *.root \n"); + return 1; + } + + TString inFileName = argv[1]; // need to check name + int timeWindow = 100; + if( argn >= 3 ) timeWindow = atoi(argv[2]); + + printf(">>> Opening input %s \n", inFileName.Data()); + TFile * inFile = new TFile(inFileName, "READ"); + if( inFile->IsOpen() == false ) { + printf("!!!! cannot open file %s \n", inFileName.Data()); + return 0; + } + + TTree * tree = (TTree *) inFile->Get("tree"); + + Long64_t evID; + UShort_t detID; + UShort_t energy; + ULong64_t energy_t; + + TBranch *b_data_ID; //! + TBranch *b_ID; //! + TBranch *b_energy; //! + TBranch *b_energy_timestamp; //! + + tree->SetBranchAddress("evID", &evID, &b_data_ID); + tree->SetBranchAddress("id", &detID, &b_ID); + tree->SetBranchAddress("e", &energy, &b_energy); + tree->SetBranchAddress("e_t", &energy_t, &b_energy_timestamp); + + Long64_t totnumEntry = tree->GetEntries(); + + printf( "total Entry : %lld \n", totnumEntry); + + printf(">>> Buidling Index using the timestamp\n"); + tree->BuildIndex("e_t"); + TTreeIndex *in = (TTreeIndex*) tree->GetTreeIndex(); + Long64_t * index = in->GetIndex(); + + ULong64_t time0; //time-0 for each event + int timeDiff; + + + TString outFileName = inFileName; + outFileName.Remove(inFileName.First("_raw")); + outFileName.Append(".root"); + if( argn >=4 ) outFileName = argv[3]; + + printf(">>> out File name : %s\n", outFileName.Data()); + + printf(">>> Create output tree\n"); + TFile * saveFile = new TFile(outFileName, "recreate"); + saveFile->cd(); + TTree * newtree = new TTree("tree", "tree"); + + newtree->Branch("evID", &eventID, "event_ID/l"); + newtree->Branch("e", e, Form("e[%d]/D", NCRYSTAL)); + newtree->Branch("e_t", e_t, Form("e_timestamp[%d]/l", NCRYSTAL)); + //newtree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCRYSTAL)); + //newtree->Branch("hit", hit, Form("hit[%d]/s", NCRYSTAL)); + + newtree->Branch("bgo", bgo, Form("BGO_e[%d]/D", NBGO)); + newtree->Branch("bgo_t", bgo_t, Form("BGO_timestamp[%d]/l", NBGO)); + + newtree->Branch("other", other, Form("other_e[%d]/D", NOTHER)); + + newtree->Branch("multi", &multi, "multiplicity_crystal/I"); + + ClearTreeData(); + + printf("================== Start processing....\n"); + Float_t Frac = 0.1; ///Progress bar + TStopwatch StpWatch; + StpWatch.Start(); + eventID = 0; + + for( Long64_t entry = 0; entry < totnumEntry; entry++){ + + + /*********** Progress Bar ******************************************/ + if (entry>totnumEntry*Frac-1) { + TString msg; msg.Form("%llu", totnumEntry/1000); + int len = msg.Sizeof(); + printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n", + Frac*100, len, entry/1000,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); + StpWatch.Start(kFALSE); + Frac+=0.1; + } + + entry = index[entry]; + + b_ID->GetEntry(entry); + b_energy->GetEntry(entry); + b_energy_timestamp->GetEntry(entry); + + if( time0 == 0) time0 = energy_t; + timeDiff = (int) (energy_t - time0); + + if( timeDiff < timeWindow ) { + + if ( detID < NCRYSTAL ){ + e[detID] = energy; + e_t[detID] = energy_t; + multi++; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = energy; + bgo_t[detID-100] = energy_t; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = energy; + } + + //printf("%d | %3d %6d %10llu, %3d\n", multi, detID, energy, energy_t, timeDiff); + + }else{ + //---- end of event + eventID ++; + + saveFile->cd(); + newtree->Fill(); + + ClearTreeData(); + + /// fill 1st data of an event + time0 = energy_t; + timeDiff = 0; + + if ( detID < NCRYSTAL ){ + e[detID] = energy; + e_t[detID] = energy_t; + multi = 1; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = energy; + bgo_t[detID-100] = energy_t; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = energy; + } + + } + + } + + printf("============================== finished.\n"); + + saveFile->cd(); + newtree->Write(); + saveFile->Close(); + + printf(" total number of event Built : %d \n", eventID); + +} diff --git a/armory/obsolete/MergeEVT.cpp b/armory/obsolete/MergeEVT.cpp new file mode 100644 index 0000000..ae58511 --- /dev/null +++ b/armory/obsolete/MergeEVT.cpp @@ -0,0 +1,132 @@ +#include +#include +#include +#include +#include +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TMath.h" +#include "TBenchmark.h" +#include + +#define MAX_CRATES 2 +#define MAX_BOARDS_PER_CRATE 13 +#define MAX_CHANNELS_PER_BOARD 16 +#define BOARD_START 2 + +#include "../mapping.h" +#include "../armory/DataBlock.h" +#include "../armory/evtReader.h" + +//############################################# +// main +//############################################# +int main(int argn, char **argv) { + + printf("=====================================\n"); + printf("=== evt --> _raw.root ===\n"); + printf("=====================================\n"); + + if (argn < 3 ) { + printf("Usage :\n"); + printf("%s [outFile] [evt1] [evt2] [evt3] ..... \n", argv[0]); + printf("e.g.: \n"); + printf("%s hahaha_raw.root haha-000.evt haha-001.evt haha-002.evt\n", argv[0]); + printf("%s hahaha_raw.root `ls haha-*.evt`\n", argv[0]); + return 1; + } + + TString outFileName = argv[1]; + int nFiles = argn-2; + TString inFileName[nFiles]; + for( int i = 0; i < nFiles ; i++){ + inFileName[i] = argv[i+2]; + printf(" in file - %2d: %s\n", i, inFileName[i].Data()); + } + + printf(" out file: %s\n", outFileName.Data()); + + + evtReader * evt = new evtReader(); + DataBlock * data = evt->data; + short detID; + printf("====================================\n"); + + //====== ROOT file + TFile * outFile = new TFile(outFileName, "recreate"); + TTree * tree = new TTree("tree", "tree"); + + tree->Branch("evID", &data->eventID, "data_ID/L"); + tree->Branch("detID", &detID, "detID/s"); + tree->Branch("e", &data->energy, "crystal_energy/s"); + tree->Branch("e_t", &data->time, "crystal_timestamp/l"); + tree->Branch("p", &data->pileup, "pileup/O"); + tree->Branch("trace_length", &data->trace_length, "trace_length/s"); + tree->Branch("trace", data->trace, "trace[trace_length]/s"); + + TBenchmark gClock; + gClock.Reset(); + gClock.Start("timer"); + + //========================================= + //========================================= + //========================================= + //========================================= + for( int i = 0; i < nFiles; i++){ + + evt->OpenFile(inFileName[i]); + if( evt->IsOpen() == false ) continue; + + Long64_t measureCount = 0; + printf("\033[1;31mProcessing file: %s\033[0m\n", inFileName[i].Data()); + TBenchmark clock2; + clock2.Reset(); + clock2.Start("timer"); + + evt->ScanNumberOfBlock(); + + //=============== Read File + while( evt->IsEndOfFile() == false ){ + + evt->ReadBlock(); + //evt->PrintStatus(10000); + + int id = data->crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data->slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data->ch; + detID = mapping[id]; + + //cern fill tree + outFile->cd(); + tree->Fill(); + } + + clock2.Stop("timer"); + double time = clock2.GetRealTime("timer"); + float tempf = (float)evt->GetFilePos()/(1024.*1024.*1024.); + printf(" measurements: \x1B[32m%lld \x1B[0m | %.3f GB\n", evt->GetBlockID(), tempf); + printf(" Time used:%3.0f min %5.2f sec\n", TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + printf(" Root file size so far: %.4f GB\n", outFile->GetSize()/1024./1024./1024.); + + } + + gClock.Stop("timer"); + double time = gClock.GetRealTime("timer"); + gClock.Start("timer"); + float tempf = (float)evt->GetFilePos()/(1024.*1024.*1024.); + printf("Total measurements: \x1B[32m%lld \x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\r", + evt->GetBlockID()+1, (100*evt->GetFilePos()/evt->GetFileSize()), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + + + //cern save root + outFile->cd(); + double totRootSize = outFile->GetSize()/1024./1024./1024.; + tree->Write(); + outFile->Close(); + + gClock.Stop("timer"); + time = gClock.GetRealTime("timer"); + printf("\n==================== finished.\r\n"); + printf("Total time spend : %3.0f min %5.2f sec\n", TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + printf(" File size of %s : %.3f GB \n", outFileName.Data(), totRootSize); + +} diff --git a/armory/obsolete/pxi-time-order.c b/armory/obsolete/pxi-time-order.c new file mode 100644 index 0000000..c625d6c --- /dev/null +++ b/armory/obsolete/pxi-time-order.c @@ -0,0 +1,412 @@ +/*********************************************************/ +/* PXI Time Order -- J.M. Allmond (ORNL) -- v1 Jul 2016 */ +/* -- v2 Feb 2018 */ +/* -- v3 Jun 2018 */ +/* -- v4 May 2019 */ +/* */ +/* !Time Order Events from Pixie-16 digitizers */ +/* !Max of: */ +/* !IDs = static, Evts = dynamic, data = dynamic */ +/* */ +/* gcc -o pxi-time-order pxi-time-order.c */ +/* ./pxi-time-order datafile */ +/*********************************************************/ + +///////////////////////////////////////////////////////// +//Code assumes that sequential sub events for a // +//specific channel are time ordered; therefore, // +//unmerge data into circular buffers on a per // +//channel id basis and then remerge channels in // +//time order. // +///////////////////////////////////////////////////////// + + +#include +#include +#include +#include +#include + +#define MAX_CRATES 2 +#define MAX_BOARDS_PER_CRATE 13 +#define MAX_CHANNELS_PER_BOARD 16 +#define BOARD_START 2 + +#define MAX_ID MAX_CRATES*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + +#define HEADER_LENGTH 4 //unit = words with 4 bytes per word +#define MAX_SUB_LENGTH 2016 //unit = words with 4 bytes per word ; 2016 --> 40 micro-second trace + 4 word header + 12 extra header + + +#define DEF_SUB_EVENTS 100 //number of events for each dynamic buffer level +#define M1_SUB_EVENTS 1000 //manual input for irregular / non-linear / non-geometric progression +#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 + +struct subevent{ + long long int timestamp; + int length; //unit = words with 4 bytes per word + unsigned int *data; +}; + +struct subevent *subevents[MAX_ID]; +int nevts[MAX_ID], iptr[MAX_ID]; +int maxevts[MAX_ID]; + + +int main(int argc, char **argv) { + + FILE *fpr, *fpw; + long int fprsize=0, fprsize_orig=0, fprsize_old=-1, fprpos=0; + + int online = 0; + + unsigned int subhead[HEADER_LENGTH]; + memset(subhead, 0, sizeof(subhead)); + int pause=0; + + long long int nwords=0, evts_tot_read=0, evts_tot_write=0; + + long long int time=0,time_old=0; + int length=0; + int chn=0; + int sln=0; + int crn=0; + int id=0; + + int idmax=0; + int totmem=0; + int outoforder=0; + int evts_old=0; + int evts_new=0; + + long long int timemin=0, timemin_old=0; + int min_id = -1; + + memset(nevts, 0, sizeof(nevts)); + memset(iptr, 0, sizeof(iptr)); + + int i=0, j=0; + div_t e_div; + + //open input event file + if ((fpr = fopen(argv[1], "r")) == NULL) { + fprintf(stderr, "Error, cannot open input file %s\n", argv[1]); + return 1; + } + + //write time order file to current location, not location of event file + char filenameto[80]; + char *filename = strrchr(argv[1], '/'); + if (filename == NULL) strcpy(filenameto,argv[1]); + else strcpy(filenameto,filename+1); + strcat(filenameto,".to"); + if ((fpw = fopen(filenameto, "w")) == NULL) { + fprintf(stderr, "Error, cannot open output file %s\n", filenameto); + return 1; + } + + //check for lockfile, active PID, and event file for auto "online" mode detection + FILE *FPLOCK; + char lockfile[1024]; + strcpy(lockfile, getenv("HOME")); + strcat(lockfile, "/.Pixie16Lock"); + int lockpid; + + FILE *FPPATH; + char pathfile[1024]; + char line[1024]; + char onlinefile[1024]; + strcpy(pathfile, getenv("HOME")); + strcat(pathfile, "/.Pixie16Path"); + + FPLOCK = fopen(lockfile, "r"); + if (FPLOCK != NULL) { + fscanf(FPLOCK, "%d", &lockpid); + fclose(FPLOCK); + //PID from lockfile matches a running PID; run timesort in "online" mode for now + if (getpgid(lockpid) >= 0) { + FPPATH = fopen(pathfile, "r"); + if (FPPATH == NULL) { + online = 0; + }else { + 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; + } + } + } + } + if (online == 1) printf("Auto Mode: \x1B[32mOnline\x1B[0m\n"); + else printf("Auto Mode: \x1B[32mOffline\x1B[0m\n"); + + //check file size for auto "online" mode + fprpos = ftell(fpr); + fseek(fpr, 0L, SEEK_END); + 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 + + + //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; + + } + + //check memory + if (totmem > MAX_MALLOC) {printf("Error: Exceeded MAX_MALLOC"); return -1;} + + //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]==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 3.5 sec) + // large gap could be from low rate or un/replug + if ( time_old == 0 || (time - time_old)/10000000 > 35 ) time_old = time; + + //fill buffers until full (online mode will stop filling buffers after 2.5 sec lag betweeen output/input) + if ( nevts[id] < maxevts[id] && ( (time - time_old)/10000000 < 25 || online == 0 ) ) { + + j = nevts[id] + iptr[id]; + if (j >= maxevts[id]) j -= maxevts[id]; + + subevents[id][j].timestamp = time; + + 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; + } + + subevents[id][j].length = 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]; + } + + nevts[id]++; + evts_tot_read++; + + } else { + pause = 1; + break; + } + + } // end while for fill buffers + ///////// + + + ///////// + // 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; + } + ///////// + + + //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); + } //end main while + + + //cleanup + fclose(fpr); + fclose(fpw); + for (i=0; i 0) printf("\x1B[31mWarning, there are %d events out of time order\x1B[0m\n", outoforder); + if (totmem != 0) printf("\x1B[31mError: total memory not conserved\x1B[0m\n"); + + + return 0; +} + + diff --git a/armory/obsolete/xia2root.cpp b/armory/obsolete/xia2root.cpp new file mode 100644 index 0000000..20037fd --- /dev/null +++ b/armory/obsolete/xia2root.cpp @@ -0,0 +1,375 @@ +#include +#include +#include +#include +#include +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TMath.h" +#include + +#define NUMDET 64 /// number of detector +#define STARTDETID 15 + +std::vector SplitStr(std::string tempLine, std::string splitter, int shift = 0){ + + std::vector 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> LoadCorrectionParameters(TString corrFile){ + + printf("==================== load correction parameters : %s", corrFile.Data()); + std::ifstream file; + file.open(corrFile.Data()); + + std::vector> corr; + corr.clear(); + + std::vector 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 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 "<= 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 "<> 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; +} +