rtang/XIAEventBuilder
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

add evt2hist.cpp

Browse Source
This commit is contained in:
Ryan Tang 2021-12-20 17:10:44 -05:00
parent c93b2a76a0
commit 0b69b17aba
6 changed files with 339 additions and 85 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -10,6 +10,7 @@ pixie2root
scan scan
pxi-time-order pxi-time-order
evt2root evt2root
evt2hist
*.so *.so
*.d *.d

66
AnalysisLibrary.h
View File

@ -231,5 +231,71 @@ std::vector<std::vector<double>> FindMatchingPair(std::vector<double> enX, std::
} }
std::vector<std::vector<double>> LoadCorrectionParameters(TString corrFile){
printf(" load correction parameters : %s", corrFile.Data());
std::ifstream file;
file.open(corrFile.Data());
std::vector<std::vector<double>> corr;
corr.clear();
std::vector<double> 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<std::string> 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<double> temp = {0, 1};
for( int i = 0; i < NCRYSTAL; i++){
corr.push_back(temp);
}
}
return corr;
}
double ApplyCorrection(std::vector<std::vector<double>> 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 #endif

58
Analyzer.h
View File

@ -113,62 +113,4 @@ void Analyzer::SlaveTerminate(){
} }
std::vector<std::vector<double>> LoadCorrectionParameters(TString corrFile){
printf(" load correction parameters : %s", corrFile.Data());
std::ifstream file;
file.open(corrFile.Data());
std::vector<std::vector<double>> corr;
corr.clear();
std::vector<double> 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<std::string> 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<double> temp = {0, 1};
for( int i = 0; i < NCRYSTAL; i++){
corr.push_back(temp);
}
}
return corr;
}
#endif // #ifdef Analyzer_cxx #endif // #ifdef Analyzer_cxx

237
evt2hist.cpp Normal file
View File

@ -0,0 +1,237 @@
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string.h>
#include <vector>
#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<std::vector<double>> 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();
}

57
evt2root.cpp
View File

@ -83,15 +83,6 @@ int main(int argn, char **argv) {
printf("====================================\n"); 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"); FILE * inFile = fopen(inFileName, "r");
if( inFile == NULL ){ if( inFile == NULL ){
printf("Cannot read file : %s \n", inFileName.Data()); printf("Cannot read file : %s \n", inFileName.Data());
@ -100,16 +91,33 @@ int main(int argn, char **argv) {
//get file size //get file size
fseek(inFile, 0L, SEEK_END); fseek(inFile, 0L, SEEK_END);
long int fprsize = ftell(inFile); long int inFileSize = ftell(inFile);
rewind(inFile); ///back to the File begining 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 int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits.
unsigned long long nWords = 0; unsigned long long nWords = 0;
unsigned long long fpos = 0;
//=============== Read File //=============== 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); fread(header, sizeof(header), 1, inFile);
fpos += sizeof(header);
measureID ++; measureID ++;
/// see the Pixie-16 user manual, Table4-2 /// see the Pixie-16 user manual, Table4-2
@ -121,7 +129,7 @@ int main(int argn, char **argv) {
data.pileup = header[0] >> 31 ; data.pileup = header[0] >> 31 ;
data.time = ((ULong64_t)(header[2] & 0xFFFF) << 32) + header[1]; data.time = ((ULong64_t)(header[2] & 0xFFFF) << 32) + header[1];
data.cfd = header[2] >> 16 ; data.cfd = header[2] >> 16 ;
data.energy = header[3] & 0xFF; data.energy = header[3] & 0xFFFF;
data.trace_length = (header[3] >> 16) & 0x7FFF; data.trace_length = (header[3] >> 16) & 0x7FFF;
data.trace_out_of_range = header[3] >> 31; data.trace_out_of_range = header[3] >> 31;
@ -129,30 +137,27 @@ int main(int argn, char **argv) {
nWords += data.eventLength; nWords += data.eventLength;
///printf("----------------------%llu\n", nWords); //printf("----------------------nWords: %llu, fpos: %llu\n", nWords, fpos);
///for(int i = 0; i < 4; i++){ //for(int i = 0; i < 4; i++){
/// printf(" %x\n", header[i]); // printf(" %x\n", header[i]);
///} //}
///data.Print(); //data.Print();
//=== jump to next measurement //=== jump to next measurement
if( data.eventLength > 4 ){ if( data.eventLength > 4 ){
fseek(inFile, sizeof(int) * (data.eventLength-4), SEEK_CUR); 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 //event stats, print status every 10000 events
if ( measureID % 10000 == 0 ) { if ( measureID % 10000 == 0 ) {
fprpos = ftell(inFile); fprpos = ftell(inFile);
float tempf = (float)fprsize/(1024.*1024.*1024.); float tempf = (float)inFileSize/(1024.*1024.*1024.);
gClock.Stop("timer"); gClock.Stop("timer");
double time = gClock.GetRealTime("timer"); double time = gClock.GetRealTime("timer");
gClock.Start("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", 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 //cern fill tree
@ -165,9 +170,9 @@ int main(int argn, char **argv) {
gClock.Stop("timer"); gClock.Stop("timer");
double time = gClock.GetRealTime("timer"); double time = gClock.GetRealTime("timer");
gClock.Start("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", 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); fclose(inFile);

5
makefile
View File

@ -1,7 +1,7 @@
CC=g++ CC=g++
#all: xia2root xia2ev2_nopart pixie2root scan pxi-time-order #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 xia2root: xia2root.cpp
$(CC) xia2root.cpp -o xia2root `root-config --cflags --glibs` $(CC) xia2root.cpp -o xia2root `root-config --cflags --glibs`
@ -15,6 +15,9 @@ pixie2root: pixie2root.cpp
evt2root: evt2root.cpp evt2root: evt2root.cpp
$(CC) evt2root.cpp -o evt2root `root-config --cflags --glibs` $(CC) evt2root.cpp -o evt2root `root-config --cflags --glibs`
evt2hist: evt2hist.cpp
$(CC) evt2hist.cpp -o evt2hist `root-config --cflags --glibs`
scan: scan.c scan: scan.c
$(CC) scan.c -o scan $(CC) scan.c -o scan