use evtReader in to2root.cpp

This commit is contained in:
Ryan Tang 2022-01-18 14:23:47 -05:00
parent c9de733c88
commit 3f8a49d83f
4 changed files with 55 additions and 167 deletions

View File

@ -62,14 +62,23 @@ public:
id = 0; id = 0;
detID = -1; detID = -1;
eventID = 0; eventID = 0;
ClearQDC();
ClearTrace();
}
void ClearQDC(){
trailing = 0; trailing = 0;
leading = 0; leading = 0;
gap = 0; gap = 0;
baseline = 0; baseline = 0;
for( int i = 0; i < 8; i++) QDCsum[i] = -1; for( int i = 0; i < 8; i++) QDCsum[i] = -1;
}
void ClearTrace(){
for( int i = 0; i < 1024; i++) trace[i] = 0; for( int i = 0; i < 1024; i++) trace[i] = 0;
} }
void Print(){ void Print(){
printf("============== eventID : %llu\n", eventID); printf("============== eventID : %llu\n", eventID);
printf("Crate: %d, Slot: %d, Ch: %d | id: %d = detID : %d \n", crate, slot, ch, id, detID); printf("Crate: %d, Slot: %d, Ch: %d | id: %d = detID : %d \n", crate, slot, ch, id, detID);

View File

@ -133,6 +133,8 @@ public:
data->id = data->crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data->slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data->ch; 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]; data->detID = mapping[data->id];
data->ClearQDC();
///======== read QDCsum ///======== read QDCsum
if( data->headerLength >= 4 ){ if( data->headerLength >= 4 ){
unsigned int extraHeader[data->headerLength-4]; unsigned int extraHeader[data->headerLength-4];
@ -150,6 +152,12 @@ public:
data->QDCsum[i] = extraHeader[i+startID]; data->QDCsum[i] = extraHeader[i+startID];
} }
} }
}else{
for( int i = 0 ; i < 8; i++){ data->QDCsum[i] = 0;}
data->trailing = 0;
data->leading = 0;
data->gap = 0;
data->baseline = 0;
} }
///====== read trace ///====== read trace
if( data->eventLength > data->headerLength ){ if( data->eventLength > data->headerLength ){

View File

@ -10,7 +10,7 @@ xia2root: ../armory/xia2root.cpp
# $(CC) armory/xia2ev2_nopart.cpp -o xia2ev2_nopart # $(CC) armory/xia2ev2_nopart.cpp -o xia2ev2_nopart
#this is for eventbuild #this is for eventbuild
to2root: ../armory/to2root.cpp to2root: ../armory/to2root.cpp ../armory/DataBlock.h ../armory/evtReader.h ../mapping.h
$(CC) ../armory/to2root.cpp -o to2root `root-config --cflags --glibs` $(CC) ../armory/to2root.cpp -o to2root `root-config --cflags --glibs`
#this is for online root #this is for online root

View File

@ -33,30 +33,7 @@
#include "../mapping.h" #include "../mapping.h"
///////////////////// #include "../armory/evtReader.h"
// RAW EVENT TYPES //
/////////////////////
struct measurement
{
int chn;
int sln;
int crn;
int id;
int hlen;
int elen;
int trlen; //number of samples
int trwlen; //number of words (two samples per word)
int fcode; //pileup flag
long long int time;
int ctime;
int ctimef;
int e;
int extra;
short tr[4096];
int esum[4];
int qsum[8];
};
struct measurement data = {0};
unsigned long long int dataCount=0; unsigned long long int dataCount=0;
unsigned long long int pileUpCount=0; unsigned long long int pileUpCount=0;
@ -67,12 +44,6 @@ unsigned long long int evtCount=0;
/////////////////////////////////// ///////////////////////////////////
int main(int argc, char **argv) { int main(int argc, char **argv) {
float tempf=0;
//temp buffer for each sub event
unsigned int sub[MAX_SUB_LENGTH];
memset(sub, 0, sizeof(sub));
printf("=====================================\n"); printf("=====================================\n");
printf("=== evt.to --> root ===\n"); printf("=== evt.to --> root ===\n");
printf("=====================================\n"); printf("=====================================\n");
@ -103,6 +74,10 @@ int main(int argc, char **argv) {
outFileName.Remove(0, pos+1); outFileName.Remove(0, pos+1);
if( argc >= 5) outFileName = argv[4]; if( argc >= 5) outFileName = argv[4];
evtReader * evt = new evtReader();
evt->OpenFile(inFileName);
DataBlock * data = evt->data;
printf(" in file : %s \n", inFileName.Data()); printf(" in file : %s \n", inFileName.Data());
printf(" our file : \033[1;31m%s\033[m\n", outFileName.Data()); printf(" our file : \033[1;31m%s\033[m\n", outFileName.Data());
@ -112,6 +87,7 @@ int main(int argc, char **argv) {
printf("------------------------ Event building time window : %d tics = %d nsec \n", timeWindow, timeWindow*10); printf("------------------------ Event building time window : %d tics = %d nsec \n", timeWindow, timeWindow*10);
TFile * outRootFile = new TFile(outFileName, "recreate"); TFile * outRootFile = new TFile(outFileName, "recreate");
outRootFile->cd(); outRootFile->cd();
TTree * tree = new TTree("tree", "tree"); TTree * tree = new TTree("tree", "tree");
@ -125,7 +101,6 @@ int main(int argc, char **argv) {
Int_t multiCry = 0 ; /// this is total multiplicity for all crystal Int_t multiCry = 0 ; /// this is total multiplicity for all crystal
//unsigned short pileup[MAXMULTI]; //unsigned short pileup[MAXMULTI];
//unsigned short hit[MAXMULTI]; // number of hit in an event
tree->Branch("evID", &evID, "event_ID/l"); tree->Branch("evID", &evID, "event_ID/l");
tree->Branch("multi", &multi, "multi/I"); tree->Branch("multi", &multi, "multi/I");
@ -135,185 +110,81 @@ int main(int argc, char **argv) {
tree->Branch("qdc", qdc, "qdc[multi][8]/I"); tree->Branch("qdc", qdc, "qdc[multi][8]/I");
tree->Branch("multiCry", &multiCry, "multiplicity_crystal/I"); tree->Branch("multiCry", &multiCry, "multiplicity_crystal/I");
//open list-mode data file from PXI digitizer
FILE *fpr = fopen(argv[1], "r");
long int fprsize,fprpos;
if ( fpr == NULL) {
fprintf(stderr, "Error, cannot open input file %s\n", argv[2]);
return 1;
}
//get file size
fseek(fpr, 0L, SEEK_END);
fprsize = ftell(fpr);
rewind(fpr);
TBenchmark gClock;
gClock.Reset();
gClock.Start("timer");
///////////////////// /////////////////////
// MAIN WHILE LOOP // // MAIN WHILE LOOP //
///////////////////// /////////////////////
while (!feof(fpr)) { //main while loop while ( evt->IsEndOfFile() == false ) { //main while loop
/////////////////////////////////
// UNPACK DATA AND EVENT BUILD //
/////////////////////////////////
long long int etime = -1; long long int etime = -1;
long long int tdif = -1; long long int tdif = -1;
while (1) { //get subevents and event build for one "event" while (1) { //get subevents and event build for one "event"
if (fread(sub, sizeof(int)*HEADER_LENGTH, 1, fpr) != 1) break; if ( evt->ReadBlock() == -1) break;
fprpos = ftell(fpr); if ( evt->GetFilePos() < evt->GetFileSize() * ( 1. - frac/100.) ) continue;
if ( fprpos < fprsize * ( 1. - frac/100.) ) {
//printf("%ld / %ld \n", fprpos, fprsize);
data.elen = (sub[0] & 0x7FFE0000) >> 17;
fseek(fpr, sizeof(int)*(data.elen - HEADER_LENGTH), SEEK_CUR);
continue;
}
data.chn = sub[0] & 0xF; /// channel in digitizer
data.sln = (sub[0] & 0xF0) >> 4; /// digitizer ID
data.crn = (sub[0] & 0xF00) >> 8; /// crate
data.id = data.crn*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data.sln-BOARD_START)*MAX_CHANNELS_PER_BOARD + data.chn;
data.hlen = (sub[0] & 0x1F000) >> 12;
data.elen = (sub[0] & 0x7FFE0000) >> 17;
data.fcode = (sub[0] & 0x80000000) >> 31;
data.time = ( (long long int)(sub[2] & 0xFFFF) << 32) + sub[1];
data.ctime = (sub[2] & 0x7FFF0000) >> 16;
data.ctimef = (sub[2] & 0x80000000) >> 31;
data.e = (sub[3] & 0xFFFF);
data.trlen = (sub[3] & 0x7FFF0000) >> 16;
data.trwlen = data.trlen / 2;
data.extra = (sub[3] & 0x80000000) >> 31;
tempf = (float)data.e/RAWE_REBIN_FACTOR;// + RAND;
data.e = (int)tempf;
if( data.hlen > 4 ){
unsigned int extraHeader[data.hlen - 4 ];
fread(extraHeader, sizeof(extraHeader), 1, fpr);
if( data.hlen == 8 || data.hlen == 16){
for( int i = 0 ; i < 4; i ++) data.esum[i] = extraHeader[i];
}
if( data.hlen == 12 || data.hlen == 16){
for( int i = 0; i < 8; i++){
int startID = 0;
if( data.hlen > 12) startID = 4; ///the 1st 4 words
data.qsum[i] = extraHeader[i+startID];
}
}
}
if( data.hlen < data.elen ){
unsigned int traceBlock[data.trlen / 2];
fread(traceBlock, sizeof(traceBlock), 1, fpr);
for( int i = 0; i < data.trlen/2 ; i++){
data.tr[2*i+0] = traceBlock[i] & 0xFFFF ;
data.tr[2*i+1] = (traceBlock[i] >> 16 ) & 0xFFFF ;
}
if( data.hlen < 12 ){
for( int i = 0 ; i < 8; i++) data.qsum[i] = 0;
for( int i = 0; i < data.trlen ; i++){
if( 0 <= i && i < 31 ) data.qsum[0] += data.tr[i];
if( 31 <= i && i < 60 ) data.qsum[1] += data.tr[i];
if( 60 <= i && i < 75 ) data.qsum[2] += data.tr[i];
if( 75 <= i && i < 95 ) data.qsum[3] += data.tr[i];
if( 95 <= i && i < 105 ) data.qsum[4] += data.tr[i];
if( 105 <= i && i < 160 ) data.qsum[5] += data.tr[i];
if( 160 <= i && i < 175 ) data.qsum[6] += data.tr[i];
if( 175 <= i && i < 200 ) data.qsum[7] += data.tr[i];
}
}
}
//Set reference time for event building //Set reference time for event building
if (etime == -1) { if (etime == -1) {
etime = data.time; etime = data->time;
tdif = 0; tdif = 0;
multi = 0; multi = 0;
}else { }else {
tdif = data.time - etime; tdif = data->time - etime;
} }
//Check for end of event, rewind, and break out of while loop //Check for end of event, rewind, and break out of while loop
if (tdif > timeWindow) { if (tdif > timeWindow) {
fseek(fpr, -sizeof(int)*data.elen, SEEK_CUR); //fwrite/fread is buffered by system ; storing this in local buffer is no faster!
etime = data->time;
tdif = 0;
multi = 0;
multiCry = 0;
id[multi] = data->detID;
e[multi] = data->energy;
e_t[multi] = data->time;
for( int i = 0; i < 8; i++) qdc[multi][i] = data->QDCsum[i];
multi++ ;
if( data->detID < 100 ) multiCry ++;
break; break;
}else{ }else{
//if within time window, fill array; //if within time window, fill array;
int detID = mapping[data.id]; id[multi] = data->detID;
id[multi] = detID; e[multi] = data->energy;
e[multi] = data.e; e_t[multi] = data->time;
e_t[multi] = data.time; for( int i = 0; i < 8; i++) qdc[multi][i] = data->QDCsum[i];
for( int i = 0; i < 8; i++) qdc[multi][i] = data.qsum[i];
multi++ ; multi++ ;
if( detID < 100 ) multiCry ++; if( data->detID < 100 ) multiCry ++;
} }
// total pileups // total pileups
if (data.fcode==1) { if (data->pileup == 1) {
pileUpCount++; pileUpCount++;
} }
//more data than just the header; read entire sub event, first rewind, then read data.elen
//fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR);
//if (fread(sub, sizeof(int)*dataBlock[sevtmult].elen, 1, fpr) != 1) break;
//if (fread(sub, sizeof(int)*data.elen, 1, fpr) != 1) break;
} //end while loop for unpacking sub events and event building for one "event" } //end while loop for unpacking sub events and event building for one "event"
if (multi==0) break; //end main WHILE LOOP when out of events if (multi==0) break; //end main WHILE LOOP when out of events
dataCount += multi; dataCount += multi;
evID ++; evID ++;
///////////////////////////////////// evt->PrintStatus(10000);
// END UNPACK DATA AND EVENT BUILD //
/////////////////////////////////////
//event stats, print status every 10000 events
if ( evID % 10000 == 0 ) {
fprpos = ftell(fpr);
tempf = (float)fprsize/(1024.*1024.*1024.);
gClock.Stop("timer");
double time = gClock.GetRealTime("timer");
gClock.Start("timer");
printf("Total dataBlock: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f <mult>)\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[3A\r",
dataCount, (int)((100*pileUpCount)/dataCount), evID+1, (float)dataCount/((float)evID+1), (100*fprpos/fprsize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
}
outRootFile->cd(); outRootFile->cd();
tree->Fill(); tree->Fill();
} // end main while loop } // end main while loop
///////////////////////// /////////////////////////
// END MAIN WHILE LOOP // // END MAIN WHILE LOOP //
///////////////////////// /////////////////////////
fprpos = ftell(fpr); evt->PrintStatus(1);
tempf = (float)fprsize/(1024.*1024.*1024.);
printf("Total SubEvents: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f <mult>)\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\n\033[3A\n",
dataCount, (int)((100*pileUpCount)/dataCount), evID+1, (float)dataCount/((float)evID+1), (100*fprpos/fprsize), tempf);
outRootFile->cd(); outRootFile->cd();
tree->Write(); tree->Write();
outRootFile->Close(); outRootFile->Close();
fclose(fpr); printf("\n\n\n==================== finished.\r\n");
gClock.Stop("timer");
double time = gClock.GetRealTime("timer");
printf("\n\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.);
return 0; return 0;
} }