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

change evt2hist and MergeEVT for evtReader

Browse Source
This commit is contained in:
Ryan Tang 2022-01-14 18:30:27 -05:00
parent 888994ae1c
commit c9de733c88
4 changed files with 145 additions and 272 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

92
armory/MergeEVT.cpp
View File

@ -17,6 +17,7 @@
#include "../mapping.h"
#include "../armory/DataBlock.h"
#include "../armory/evtReader.h"
//#############################################
// main
@ -46,8 +47,9 @@ int main(int argn, char **argv) {
printf(" out file: %s\n", outFileName.Data());
Long64_t measureID = -1;
DataBlock data;
evtReader * evt = new evtReader();
DataBlock * data = evt->data;
printf("====================================\n");
@ -55,29 +57,25 @@ int main(int argn, char **argv) {
TFile * outFile = new TFile(outFileName, "recreate");
TTree * tree = new TTree("tree", "tree");
tree->Branch("evID", &measureID, "data_ID/L");
tree->Branch("id", &data.id, "ID/s");
tree->Branch("e", &data.energy, "crystal_energy/s");
tree->Branch("e_t", &data.time, "crystal_timestamp/l");
tree->Branch("evID", &data->eventID, "data_ID/L");
tree->Branch("id", &data->id, "ID/s");
tree->Branch("e", &data->energy, "crystal_energy/s");
tree->Branch("e_t", &data->time, "crystal_timestamp/l");
long int inFileSize = 0;
long int inFilePos = 0;
TBenchmark gClock;
gClock.Reset();
gClock.Start("timer");
//=========================================
//=========================================
//=========================================
//=========================================
for( int i = 0; i < nFiles; i++){
FILE * inFile = fopen(inFileName[i], "r");
if( inFile == NULL ){
printf("Cannot read file : %s \n", inFileName[i].Data());
continue;
}
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());
@ -85,74 +83,22 @@ int main(int argn, char **argv) {
clock2.Reset();
clock2.Start("timer");
//get file size
fseek(inFile, 0L, SEEK_END);
inFileSize = ftell(inFile);
rewind(inFile); ///back to the File begining
inFilePos = 0;
unsigned int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits.
unsigned long long nWords = 0;
//ULong64_t timeLast = 0;
//=============== Read File
/// while ( ! feof(inFile) ){
while( inFilePos < inFileSize || feof(inFile) ){
fread(header, sizeof(header), 1, inFile);
inFilePos = ftell(inFile);
measureID ++;
measureCount++;
while( evt->IsEndOfFile() == false ){
/// 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;
data.id = mapping[data.id];
data.energy = data.energy / 2. ; // factor 2 is the rawe_rebin_factor;
nWords += data.eventLength;
//=== jump to next measurement
if( data.eventLength > 4 ){
fseek(inFile, sizeof(int) * (data.eventLength-4), SEEK_CUR);
inFilePos = ftell(inFile);
}
//event stats, print status every 10000 events
if ( measureID % 10000 == 0 ) {
float tempf = (float)inFileSize/(1024.*1024.*1024.);
gClock.Stop("timer");
double time = gClock.GetRealTime("timer");
gClock.Start("timer");
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\033[A\r",
measureID +1 , (100*inFilePos/inFileSize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
}
evt->ReadBlock();
evt->PrintStatus(10000);
//cern fill tree
outFile->cd();
tree->Fill();
}
inFilePos = ftell(inFile);
fclose(inFile);
clock2.Stop("timer");
double time = clock2.GetRealTime("timer");
float tempf = (float)inFileSize/(1024.*1024.*1024.);
printf(" measurements: \x1B[32m%lld \x1B[0m | %.3f GB\n", measureCount, tempf);
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.);
@ -161,9 +107,9 @@ int main(int argn, char **argv) {
gClock.Stop("timer");
double time = gClock.GetRealTime("timer");
gClock.Start("timer");
float tempf = (float)inFileSize/(1024.*1024.*1024.);
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",
measureID+1, (100*inFilePos/inFileSize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
evt->GetBlockID()+1, (100*evt->GetFilePos()/evt->GetFileSize()), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
//cern save root

240
armory/evt2hist.cpp
View File

@ -16,6 +16,7 @@
#include "TLatex.h"
#include "TBenchmark.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TApplication.h"
#include "TCanvas.h"
#include "TClonesArray.h"
@ -24,28 +25,14 @@
#include "../armory/AnalysisLibrary.h"
#include "../armory/DataBlock.h"
#include "../armory/evtReader.h"
#define MAX_CRATES 2
#define MAX_BOARDS_PER_CRATE 13
#define MAX_CHANNELS_PER_BOARD 16
#define BOARD_START 2
#define sleepTime 5000 ///sleep for 5 sec
//#############################TODO
// 1) multiple file
// 2) Change to GUI
// 4) eventBuilding
DataBlock dataList[1000];
int dataCollected = 0;
void BuildEvent(){
///==== sort timestamp
///==== build events
///==== output to g-g hist
}
// 3) last 10% data
//#############################################
// main
@ -99,22 +86,17 @@ int main(int argn, char **argv) {
int maxDataDisplay = 0;
if( argn >= 7 ) maxDataDisplay = atoi(argv[6]);
long int inFilePos;
TBenchmark gClock;
gClock.Reset();
gClock.Start("timer");
ULong64_t measureID = -1;
DataBlock data;
printf("====================================\n");
FILE * inFile = fopen(inFileName, "r");
if( inFile == NULL ){
printf("Cannot read file : %s \n", inFileName.Data());
return -404;
}
evtReader * evt = new evtReader();
evt->OpenFile(inFileName);
if( evt->IsOpen() == false ) return -404;
DataBlock * data = evt->data;
printf(" in file: \033[1;31m%s\033[m\n", inFileName.Data());
printf(" Gamma energy correction file : %s\n", corrFile == "" ? "Not provided." : corrFile.Data());
@ -123,6 +105,8 @@ int main(int argn, char **argv) {
if( rootFileName != "" ) printf(" Save root to %s\n", rootFileName.Data());
printf("--------------------------------\n");
//================ ROOT tree
TFile * fFile = NULL;
TTree * tree = NULL;
if( rootFileName != "" ){
@ -130,23 +114,17 @@ int main(int argn, char **argv) {
tree = new TTree("tree", "tree");
tree->Branch("headerLenght", &data.headerLength, "HeaderLength/s");
tree->Branch("detID", &data.detID, "detID/s");
tree->Branch("id", &data.id, "id/s");
tree->Branch("e", &data.energy, "energy/s");
tree->Branch("e_t", &data.time, "timestamp/l");
tree->Branch("p", &data.pileup, "pileup/O");
tree->Branch("qdc", data.QDCsum, "QDC_sum[8]/I");
tree->Branch("trace_length", &data.trace_length, "trace_length/s");
tree->Branch("trace", data.trace, "trace[trace_length]/s");
tree->Branch("headerLenght", &data->headerLength, "HeaderLength/s");
tree->Branch("detID", &data->detID, "detID/s");
tree->Branch("id", &data->id, "id/s");
tree->Branch("e", &data->energy, "energy/s");
tree->Branch("e_t", &data->time, "timestamp/l");
tree->Branch("p", &data->pileup, "pileup/O");
tree->Branch("qdc", data->QDCsum, "QDC_sum[8]/I");
tree->Branch("trace_length", &data->trace_length, "trace_length/s");
tree->Branch("trace", data->trace, "trace[trace_length]/s");
}
//================ get file size
fseek(inFile, 0L, SEEK_END);
long int inFileSize = ftell(inFile);
rewind(inFile); ///back to the File begining
unsigned long long fpos = 0;
//================ Historgrams
TH1F * he[NCRYSTAL];
for( int i = 0 ; i < NCRYSTAL; i++){
@ -159,6 +137,9 @@ int main(int argn, char **argv) {
}
}
int GAGGID = 209;
TH2F * hPID = new TH2F(Form("hPID%d", GAGGID), Form("GAGG - %d; tail; peak", GAGGID), 400, -10, 600, 400, -50, 1000);
TGraph * gTrace = new TGraph();
TLatex text;
text.SetNDC();
@ -180,116 +161,69 @@ int main(int argn, char **argv) {
///TCanvas * cTrace = new TCanvas("cTrace", "Trace", 100, 100, 1000, 500);
TCanvas * cPID = new TCanvas("cPID", "PID", 100, 100, 500, 500);
//=============== Read File
unsigned int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits.
int sleepCount = 0;
while ( true ){
while ( ! feof(inFile) ){
if ( fread(header, sizeof(header), 1, inFile) !=1 ) break;
measureID ++;
int status = evt->ReadBlock();
/// see the Pixie-16 user manual, Table4-2
data.eventID = measureID;
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 ) /2; // I don;t know why it has to "rebin"
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;
data.detID = mapping[data.id];
///======== read QDCsum
if( data.headerLength >= 4 ){
unsigned int extraHeader[data.headerLength-4];
fread(extraHeader, sizeof(extraHeader), 1, inFile);
if( data.headerLength == 8 || data.headerLength == 16){
data.trailing = extraHeader[0];
data.leading = extraHeader[1];
data.gap = extraHeader[2];
data.baseline = extraHeader[3];
}
if( data.headerLength == 12 || data.headerLength == 16){
for( int i = 0; i < 8; i++){
int startID = 0;
if( data.headerLength > 12) startID = 4; ///the 1st 4 words
data.QDCsum[i] = extraHeader[i+startID];
}
}
if( status == -1 ) {
break;
//printf("\n\n\nReached the end of file, wait %d sec to see any update.\n", sleepTime);
//sleep( sleepTime );
//evt->UpdateFileSize();
//sleepCount ++;
//if( sleepCount > 1 ) {
// printf("waited for %d sec. exit.\n", 4* sleepTime);
// break;
//}else{
// continue;
//}
}
///====== read trace
if( data.eventLength > data.headerLength ){
unsigned int traceBlock[data.trace_length / 2];
fread(traceBlock, sizeof(traceBlock), 1, inFile);
for( int i = 0; i < data.trace_length/2 ; i++){
data.trace[2*i+0] = traceBlock[i] & 0xFFFF ;
data.trace[2*i+1] = (traceBlock[i] >> 16 ) & 0xFFFF ;
}
sleepCount = 0;
if( evt->GetBlockID() < maxDataDisplay ) {
printf("----------------------event Length: %u, fpos: %lu byte (%ld words)\n", data->eventLength, evt->GetFilePos(), evt->GetFilePos()/4);
data->Print();
}
if( measureID < maxDataDisplay ) {
printf("----------------------event Length: %u, fpos: %llu byte (%lld words)\n", data.eventLength, fpos, fpos/4);
for(int i = 0; i < 4; i++) printf(" %x\n", header[i]);
data.Print();
}
//=== jump to next measurement. obsolete, we read the whole block
///if( data.eventLength > 4 ) {
/// if( fseek(inFile, sizeof(int) * (data.eventLength-data.headerLength), SEEK_CUR) != 0 ) break;;
///}
fpos = ftell(inFile);
//==== Fill Histogram
if( 0 <= data.detID && data.detID < 100 && data.energy > rawEnergyThreshold ){
if( 0 <= data->detID && data->detID < 100 && data->energy > rawEnergyThreshold ){
if( corrFile != ""){
///========= apply correction
int order = (int) eCorr[data.detID].size();
int order = (int) eCorr[data->detID].size();
double eCal = 0;
for( int k = 0; k < order ; k++){
eCal += eCorr[data.detID][k] * TMath::Power(data.energy, k);
eCal += eCorr[data->detID][k] * TMath::Power(data->energy, k);
}
he[data.detID]->Fill(eCal);
he[data->detID]->Fill(eCal);
}else{
he[data.detID]->Fill(data.energy);
he[data->detID]->Fill(data->energy);
}
}
///============ QDC
if( data->detID == GAGGID && (data->headerLength < data->eventLength) ){
double bg = (data->QDCsum[0] + data->QDCsum[1])/60.;
double peak = data->QDCsum[3]/20. - bg;
double tail = data->QDCsum[5]/55. - bg;
hPID->Fill( tail , peak);
}
//===== Trace
if( data.trace_length > 0 ) {
gTrace->Clear();
gTrace->Set(data.trace_length);
gTrace->SetTitle(Form("eventID : %llu, detID: %d", data.eventID, data.detID));
if( data.headerLength == 4 ) {
for( int i = 0; i < 8; i++ ) data.QDCsum[i] = 0;
}
for( int i = 0; i < data.trace_length; i++){
gTrace->SetPoint(i, i, data.trace[i]);
///if the header don't have ADC, make one
if( data.headerLength < 12 ) {
if( 0 <= i && i < 31 ) data.QDCsum[0] += data.trace[i];
if( 31 <= i && i < 60 ) data.QDCsum[1] += data.trace[i];
if( 60 <= i && i < 75 ) data.QDCsum[2] += data.trace[i];
if( 75 <= i && i < 95 ) data.QDCsum[3] += data.trace[i];
if( 95 <= i && i < 105 ) data.QDCsum[4] += data.trace[i];
if( 105 <= i && i < 160 ) data.QDCsum[5] += data.trace[i];
if( 160 <= i && i < 175 ) data.QDCsum[6] += data.trace[i];
if( 175 <= i && i < 200 ) data.QDCsum[7] += data.trace[i];
}
}
}
///if( data->trace_length > 0 ) {
/// gTrace->Clear();
/// gTrace->Set(data->trace_length);
/// gTrace->SetTitle(Form("eventID : %llu, detID: %d", data->eventID, data->detID));
///
/// for( int i = 0; i < data->trace_length; i++) gTrace->SetPoint(i, i, data->trace[i]);
///}
if( rootFileName != "" ){
fFile->cd();
@ -297,19 +231,7 @@ int main(int argn, char **argv) {
}
//==== event stats, print status every 10000 events
if ( measureID % 10000 == 0 ) {
/// update file size, slow down?
fseek(inFile, 0L, SEEK_END);
inFileSize = ftell(inFile);
fseek(inFile, fpos, SEEK_SET);
inFilePos = ftell(inFile);
gClock.Stop("timer");
double time = gClock.GetRealTime("timer");
gClock.Start("timer");
printf("Total measurements: \x1B[32m%llu \x1B[0m\nReading Pos: \x1B[32m %.3f/%.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\033[A\r",
measureID, inFilePos/(1024.*1024.*1024.), inFileSize/1024./1024./1024, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
}
evt->PrintStatus(10000);
//==== Plot Canvas
gClock.Stop("timer");
@ -317,7 +239,7 @@ int main(int argn, char **argv) {
gClock.Start("timer");
if( time % 1000 == 0 || time < 10){
//==== for clover
///==== for clover
for( int i = 0; i < NCLOVER; i++){
double maxY = 0;
double y = 0;
@ -339,18 +261,24 @@ int main(int argn, char **argv) {
canvas->Modified();
canvas->Update();
//==== for trace
///if( data.trace_length > 0 ){
///==== for trace
///if( data->trace_length > 0 ){
/// cTrace->cd();
/// gTrace->Draw("AL");
///
/// for( int i = 0; i < 8; i++){
/// text.DrawLatex(0.2, 0.8-0.05*i, Form("%d", data.QDCsum[i]));
/// text.DrawLatex(0.2, 0.8-0.05*i, Form("%d", data->QDCsum[i]));
/// }
/// cTrace->Modified();
/// cTrace->Update();
///}
///=== for GAGG PID
cPID->cd();
cPID->SetLogz();
hPID->Draw("colz");
cPID->Modified();
cPID->Update();
gSystem->ProcessEvents();
}
@ -380,15 +308,7 @@ int main(int argn, char **argv) {
gSystem->ProcessEvents();
inFilePos = ftell(inFile);
gClock.Stop("timer");
double time = gClock.GetRealTime("timer");
gClock.Start("timer");
printf("Total measurements: \x1B[32m%llu \x1B[0m\nReading Pos: \x1B[32m %.3f/%.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\033[A\r",
measureID, inFilePos/(1024.*1024.*1024.), inFileSize/1024./1024./1024, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
fclose(inFile);
evt->PrintStatus(1);
printf("\n\n\n============= reached end of file\n");

79
armory/evtReader.h
View File

@ -27,53 +27,34 @@ class evtReader{
public:
DataBlock * data;
Long64_t blockID;
bool isOpen;
private:
FILE * inFile;
long int inFileSize;
long int inFilePos;
bool endOfFile;
bool isOpened;
Long64_t blockID;
TBenchmark gClock;
///============================================ Methods
public:
evtReader(TString inFileName){
inFile = fopen(inFileName, "r");
if( inFile == NULL ){
printf("Cannot read file : %s \n", inFileName.Data());
inFileSize = 0;
inFilePos = 0;
blockID = -1;
endOfFile = false;
isOpen = false;
}else{
evtReader(){
inFile = 0;
data = new DataBlock();
fseek(inFile, 0L, SEEK_END);
inFileSize = ftell(inFile);
inFilePos = 0;
rewind(inFile); ///back to the File begining
inFileSize = 0;
inFilePos = 0;
data = new DataBlock();
data->Clear();
blockID = -1;
endOfFile = false;
isOpen = true;
gClock.Reset();
gClock.Start("timer");
}
blockID = -1;
endOfFile = false;
isOpened = false;
}
evtReader(TString inFileName){
OpenFile(inFileName);
}
~evtReader(){
fclose(inFile);
@ -81,6 +62,28 @@ public:
delete data;
}
void OpenFile(TString inFileName){
inFile = fopen(inFileName, "r");
if( inFile == NULL ){
printf("Cannot read file : %s \n", inFileName.Data());
}else{
fseek(inFile, 0L, SEEK_END);
inFileSize = ftell(inFile);
inFilePos = 0;
rewind(inFile); ///back to the File begining
data->Clear();
blockID = -1;
endOfFile = false;
isOpened = true;
gClock.Reset();
gClock.Start("timer");
}
}
void UpdateFileSize(){
if( inFile == NULL ) return;
fseek(inFile, 0L, SEEK_END);
@ -88,7 +91,9 @@ public:
fseek(inFile, inFilePos, SEEK_SET);
}
bool isEndOfFile() {
bool IsOpen(){ return isOpened;}
bool IsEndOfFile() {
int haha = feof(inFile);
return haha > 0 ? true: false;
}
@ -96,6 +101,8 @@ public:
long int GetFilePos(){return inFilePos;}
long int GetFileSize(){return inFileSize;}
Long64_t GetBlockID(){ return blockID;}
int ReadBlock(){
if( feof(inFile) ) return -1;
@ -155,7 +162,8 @@ public:
}
///make QDC by trace
if( data->headerLength == 8 ) {
if( data->headerLength == 4 || data->headerLength == 8 ) {
for( int i = 0; i < 8; i++){ data->QDCsum[i] = 0;}
for( int i = 0; i < data->trace_length; i++){
if( 0 <= i && i < 31 ) data->QDCsum[0] += data->trace[i];
if( 31 <= i && i < 60 ) data->QDCsum[1] += data->trace[i];
@ -180,7 +188,6 @@ public:
///==== event stats, print status every 10000 events
if ( blockID % id == 0 ) {
UpdateFileSize();
//inFilePos = ftell(inFile);
gClock.Stop("timer");
double time = gClock.GetRealTime("timer");
gClock.Start("timer");

6
armory/makefile
View File

@ -1,6 +1,6 @@
CC=g++
all: to2root evt2hist MergeEVT ev22txt EventBuilder pxi-time-order test
all: to2root evt2hist MergeEVT ev22txt EventBuilder pxi-time-order
#this is FSU evt to root
xia2root: ../armory/xia2root.cpp
@ -14,11 +14,11 @@ to2root: ../armory/to2root.cpp
$(CC) ../armory/to2root.cpp -o to2root `root-config --cflags --glibs`
#this is for online root
MergeEVT: ../armory/MergeEVT.cpp ../armory/DataBlock.h ../mapping.h
MergeEVT: ../armory/MergeEVT.cpp ../armory/DataBlock.h ../armory/evtReader.h ../mapping.h
$(CC) ../armory/MergeEVT.cpp -o MergeEVT `root-config --cflags --glibs`
#this is for online spectrums
evt2hist: ../armory/evt2hist.cpp ../armory/DataBlock.h ../mapping.h
evt2hist: ../armory/evt2hist.cpp ../armory/DataBlock.h ../armory/evtReader.h ../mapping.h
$(CC) ../armory/evt2hist.cpp -o evt2hist `root-config --cflags --glibs`
pxi-time-order: ../armory/pxi-time-order.c