From 0b69b17aba42cffb7c37c943c13703d40833879e Mon Sep 17 00:00:00 2001 From: "Ryan@WorkStation" Date: Mon, 20 Dec 2021 17:10:44 -0500 Subject: [PATCH] add evt2hist.cpp --- .gitignore | 1 + AnalysisLibrary.h | 66 +++++++++++++ Analyzer.h | 58 ------------ evt2hist.cpp | 237 ++++++++++++++++++++++++++++++++++++++++++++++ evt2root.cpp | 57 ++++++----- makefile | 5 +- 6 files changed, 339 insertions(+), 85 deletions(-) create mode 100644 evt2hist.cpp diff --git a/.gitignore b/.gitignore index 1245ef3..b889b66 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ pixie2root scan pxi-time-order evt2root +evt2hist *.so *.d diff --git a/AnalysisLibrary.h b/AnalysisLibrary.h index baf2118..030a6f6 100644 --- a/AnalysisLibrary.h +++ b/AnalysisLibrary.h @@ -231,5 +231,71 @@ std::vector> FindMatchingPair(std::vector enX, std:: } + +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("%14.6f, ", corr[i][j]); + } + printf("%14.6f\n", corr[i][len-1]); + } + + }else{ + printf(".... fail\n"); + std::vector temp = {0, 1}; + for( int i = 0; i < NCRYSTAL; i++){ + corr.push_back(temp); + } + } + + return corr; +} + +double ApplyCorrection(std::vector> corr, int corrRow, double value){ + double eCal = 0 ; + int order = (int) corr[corrRow].size(); + for( int i = 0; i < order ; i++){ + eCal += corr[corrRow][i] * TMath::Power(value, i); + } + return eCal; +} + #endif diff --git a/Analyzer.h b/Analyzer.h index 6d3d8e3..9304948 100644 --- a/Analyzer.h +++ b/Analyzer.h @@ -113,62 +113,4 @@ void Analyzer::SlaveTerminate(){ } -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("%14.6f, ", corr[i][j]); - } - printf("%14.6f\n", corr[i][len-1]); - } - - }else{ - printf(".... fail\n"); - std::vector temp = {0, 1}; - for( int i = 0; i < NCRYSTAL; i++){ - corr.push_back(temp); - } - } - - return corr; -} - - - #endif // #ifdef Analyzer_cxx diff --git a/evt2hist.cpp b/evt2hist.cpp new file mode 100644 index 0000000..fa1b10e --- /dev/null +++ b/evt2hist.cpp @@ -0,0 +1,237 @@ +#include +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TMath.h" +#include "TBenchmark.h" +#include "TH1F.h" +#include "TApplication.h" +#include "TCanvas.h" +#include "TSystem.h" + +#include "mapping.h" +#include "AnalysisLibrary.h" + +#define MAX_CRATES 2 +#define MAX_BOARDS_PER_CRATE 13 +#define MAX_CHANNELS_PER_BOARD 16 +#define BOARD_START 2 + + +class measurment{ + +public: + UShort_t ch; + UShort_t slot; + UShort_t crate; + UShort_t headerLength; /// headerLength > 4, more data except tarce. + UShort_t eventLength; /// eventLength = headerLength + trace + Bool_t pileup; + ULong64_t time; + UShort_t cfd; + UShort_t energy; + UShort_t trace_length; + Bool_t trace_out_of_range; + + UShort_t id; + + measurment(){}; + + void Clear(){ + ch = 0; + slot = 0; + crate = 0; + eventLength = 0; + pileup = false; + time = 0; + cfd = 0; + energy = 0; + trace_length = 0; + trace_out_of_range = 0; + id = 0; + } + + void Print(){ + printf("Crate: %d, Slot: %d, Ch: %d | id: %d \n", crate, slot, ch, id); + printf("HeaderLength: %d, Event Length: %d, energy: %d, timeStamp: %llu\n", headerLength, eventLength, energy, time); + printf("trace_length: %d, pile-up:%d\n", trace_length, pileup); + } + +}; + +//############################################# +// main +//############################################### +int main(int argn, char **argv) { + + if (argn != 2 && argn != 3 ) { + printf("Usage :\n"); + printf("%s [evt File] [E corr]\n", argv[0]); + printf(" [E corr] : correction file for gamma energy \n"); + return 1; + } + + TString inFileName = argv[1]; + + TString corrFile = ""; + std::vector> corr; + if( argn >= 3 ){ + corrFile = argv[2]; + corr.clear(); + corr = LoadCorrectionParameters(corrFile); + } + + long int inFilePos; + TBenchmark gClock; + gClock.Reset(); + gClock.Start("timer"); + + ULong64_t measureID = 0; + + measurment data; + + printf("====================================\n"); + + FILE * inFile = fopen(inFileName, "r"); + if( inFile == NULL ){ + printf("Cannot read file : %s \n", inFileName.Data()); + return -404; + } + + //get file size + //fseek(inFile, 0L, SEEK_END); + //long int inFileSize = ftell(inFile); + //rewind(inFile); ///back to the File begining + + printf(" in file: %s\n", inFileName.Data()); + printf(" Gamma energy correction file : %s\n", corrFile == "" ? "Not provided." : corrFile.Data()); + printf("--------------------------------\n"); + + //================ Historgrams + TH1F * he[NCRYSTAL]; + for( int i = 0 ; i < NCRYSTAL; i++){ + he[i] = new TH1F(Form("he%02d", i), Form("e-%2d", i), 1000, 0, 6000); + switch (i % 4){ + case 0 : he[i]->SetLineColor(2); break; + case 1 : he[i]->SetLineColor(4); break; + case 2 : he[i]->SetLineColor(1); break; + case 3 : he[i]->SetLineColor(kGreen+3); break; + } + } + + //================ Set Canvas + TApplication * app = new TApplication ("app", &argn, argv); + TCanvas * canvas = new TCanvas("fCanvas", "Online Spectrum", 1800, 1400); + canvas->Divide(1, 9, 0); + canvas->SetCrosshair(1); + for( int i = 0; i < 9 ; i++){ + canvas->cd(i+1)->SetBottomMargin(0.06); + canvas->cd(i+1)->SetRightMargin(0.002); + } + + unsigned int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits. + unsigned long long nWords = 0; + unsigned long long fpos = 0; + + //=============== Read File + while ( ! feof(inFile) ){ + + fread(header, sizeof(header), 1, inFile); + measureID ++; + + /// see the Pixie-16 user manual, Table4-2 + data.ch = header[0] & 0xF ; + data.slot = (header[0] >> 4) & 0xF; + data.crate = (header[0] >> 8) & 0xF; + data.headerLength = (header[0] >> 12) & 0x1F; + data.eventLength = (header[0] >> 17) & 0x3FFF; + data.pileup = header[0] >> 31 ; + data.time = ((ULong64_t)(header[2] & 0xFFFF) << 32) + header[1]; + data.cfd = header[2] >> 16 ; + data.energy = header[3] & 0xFFFF; + data.trace_length = (header[3] >> 16) & 0x7FFF; + data.trace_out_of_range = header[3] >> 31; + + data.id = data.crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data.slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data.ch; + + nWords += data.eventLength; + + ///printf("----------------------nWords: %llu, fpos: %llu\n", nWords, fpos); + ///for(int i = 0; i < 4; i++){ + /// printf(" %x\n", header[i]); + ///} + ///data.Print(); + + int detID = map[data.id]; + if( 0 <= detID && detID < 100 ){ + if( corrFile != ""){ + double eCal = ApplyCorrection(corr, detID, data.energy); + he[detID]->Fill(eCal); + }else{ + he[detID]->Fill(data.energy); + } + } + + //=== jump to next measurement + if( data.eventLength > 4 ){ + fseek(inFile, sizeof(int) * (data.eventLength-4), SEEK_CUR); + fpos += ftell(inFile); + } + + //==== event stats, print status every 10000 events + if ( measureID % 10000 == 0 ) { + //fseek(inFile, 0L, SEEK_END); + //inFileSize = ftell(inFile); + //fseek(inFile, fpos, SEEK_SET); + + inFilePos = ftell(inFile); + float tempf = (float)inFilePos/(1024.*1024.*1024.); + gClock.Stop("timer"); + double time = gClock.GetRealTime("timer"); + gClock.Start("timer"); + printf("Total measurements: \x1B[32m%llu \x1B[0m\nData Read: \x1B[32m %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\033[A\r", + measureID, tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + } + + + gClock.Stop("timer"); + int time = TMath::Floor(gClock.GetRealTime("timer")*1000); // in millisec + gClock.Start("timer"); + if( time % 1000 == 0 ){ + for( int i = 0; i < NCRYSTAL; i++){ + canvas->cd(i/4 +1); + canvas->cd(i/4 +1)->SetLogy(); + if( i % 4 == 0 ) { + he[i]->Draw(); + }else{ + he[i]->Draw("same"); + } + } + canvas->Modified(); + canvas->Update(); + gSystem->ProcessEvents(); + } + } + + inFilePos = ftell(inFile); + gClock.Stop("timer"); + double time = gClock.GetRealTime("timer"); + gClock.Start("timer"); + float tempf = (float)inFilePos/(1024.*1024.*1024.); + printf("Total measurements: \x1B[32m%llu \x1B[0m\nData Read: \x1B[32m %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\033[A\r", + measureID, tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + + fclose(inFile); + + printf("\n============= reasched end of file\n"); + printf("Crtl+C to end program.\n"); + app->Run(); + + +} diff --git a/evt2root.cpp b/evt2root.cpp index c9d940d..d69ff54 100644 --- a/evt2root.cpp +++ b/evt2root.cpp @@ -83,15 +83,6 @@ int main(int argn, char **argv) { printf("====================================\n"); - TFile * outFile = new TFile(outFileName, "recreate"); - TTree * tree = new TTree("tree", "tree"); - - tree->Branch("evID", &measureID, "data_ID/l"); - tree->Branch("detID", &data.id, "det_ID/s"); - tree->Branch("e", &data.energy, "energy/s"); - tree->Branch("t", &data.time, "time_stamp/l"); - - FILE * inFile = fopen(inFileName, "r"); if( inFile == NULL ){ printf("Cannot read file : %s \n", inFileName.Data()); @@ -100,16 +91,33 @@ int main(int argn, char **argv) { //get file size fseek(inFile, 0L, SEEK_END); - long int fprsize = ftell(inFile); + long int inFileSize = ftell(inFile); rewind(inFile); ///back to the File begining + printf(" in file: %s\n", inFileName.Data()); + printf("out file: %s\n", outFileName.Data()); + printf("--------------------------------\n"); + + TFile * outFile = new TFile(outFileName, "recreate"); + TTree * tree = new TTree("tree", "tree"); + + tree->Branch("evID", &measureID, "data_ID/l"); + tree->Branch("detID", &data.id, "det_ID/s"); + tree->Branch("e", &data.energy, "energy/s"); + tree->Branch("t", &data.time, "time_stamp/l"); + +//=======TODO online event building + unsigned int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits. unsigned long long nWords = 0; + unsigned long long fpos = 0; //=============== Read File - while ( ! feof(inFile) ){ +// while ( ! feof(inFile) ){ + while ( fpos <= inFileSize ){ // need to check is the last data included. fread(header, sizeof(header), 1, inFile); + fpos += sizeof(header); measureID ++; /// see the Pixie-16 user manual, Table4-2 @@ -121,7 +129,7 @@ int main(int argn, char **argv) { data.pileup = header[0] >> 31 ; data.time = ((ULong64_t)(header[2] & 0xFFFF) << 32) + header[1]; data.cfd = header[2] >> 16 ; - data.energy = header[3] & 0xFF; + data.energy = header[3] & 0xFFFF; data.trace_length = (header[3] >> 16) & 0x7FFF; data.trace_out_of_range = header[3] >> 31; @@ -129,30 +137,27 @@ int main(int argn, char **argv) { nWords += data.eventLength; - ///printf("----------------------%llu\n", nWords); - ///for(int i = 0; i < 4; i++){ - /// printf(" %x\n", header[i]); - ///} - ///data.Print(); - + //printf("----------------------nWords: %llu, fpos: %llu\n", nWords, fpos); + //for(int i = 0; i < 4; i++){ + // printf(" %x\n", header[i]); + //} + //data.Print(); //=== jump to next measurement if( data.eventLength > 4 ){ fseek(inFile, sizeof(int) * (data.eventLength-4), SEEK_CUR); + fpos += sizeof(int) * (data.eventLength-4); } - - //if( nWords > 200 ) break; - - + //event stats, print status every 10000 events if ( measureID % 10000 == 0 ) { fprpos = ftell(inFile); - float tempf = (float)fprsize/(1024.*1024.*1024.); + float tempf = (float)inFileSize/(1024.*1024.*1024.); gClock.Stop("timer"); double time = gClock.GetRealTime("timer"); gClock.Start("timer"); printf("Total measurements: \x1B[32m%llu \x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\033[A\r", - measureID, (100*fprpos/fprsize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + measureID, (100*fprpos/inFileSize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); } //cern fill tree @@ -165,9 +170,9 @@ int main(int argn, char **argv) { gClock.Stop("timer"); double time = gClock.GetRealTime("timer"); gClock.Start("timer"); - float tempf = (float)fprsize/(1024.*1024.*1024.); + float tempf = (float)inFileSize/(1024.*1024.*1024.); printf("Total measurements: \x1B[32m%llu \x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\r", - measureID, (100*fprpos/fprsize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + measureID, (100*fprpos/inFileSize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); fclose(inFile); diff --git a/makefile b/makefile index 44b08d1..5f02089 100644 --- a/makefile +++ b/makefile @@ -1,7 +1,7 @@ CC=g++ #all: xia2root xia2ev2_nopart pixie2root scan pxi-time-order -all: xia2root xia2ev2_nopart pixie2root scan evt2root +all: xia2root xia2ev2_nopart pixie2root scan evt2root evt2hist xia2root: xia2root.cpp $(CC) xia2root.cpp -o xia2root `root-config --cflags --glibs` @@ -15,6 +15,9 @@ pixie2root: pixie2root.cpp evt2root: evt2root.cpp $(CC) evt2root.cpp -o evt2root `root-config --cflags --glibs` +evt2hist: evt2hist.cpp + $(CC) evt2hist.cpp -o evt2hist `root-config --cflags --glibs` + scan: scan.c $(CC) scan.c -o scan