Make EventBuilder, skip the *.to file. clean up armory

This commit is contained in:
Ryan Tang 2024-05-29 14:39:38 -04:00
parent b44e63377a
commit a727af0b48
14 changed files with 1108 additions and 2165 deletions

1
.gitignore vendored
View File

@ -15,6 +15,7 @@ MergeEVT
evt2hist
ev22txt
EventBuilder
EventBuilder2
*.so
*.d

View File

@ -14,6 +14,17 @@
"cStandard": "c17",
"cppStandard": "c++17",
"intelliSenseMode": "macos-clang-arm64"
},
{
"name": "Ryan",
"includePath": [
"${workspaceFolder}/**",
"/opt/root/**"
],
"defines": [],
"compilerPath": "/usr/bin/clang",
"cStandard": "c17",
"cppStandard": "c++17"
}
],
"version": 4

31
.vscode/settings.json vendored
View File

@ -64,6 +64,35 @@
"unordered_set": "cpp",
"valarray": "cpp",
"variant": "cpp",
"vector": "cpp"
"vector": "cpp",
"any": "cpp",
"atomic": "cpp",
"strstream": "cpp",
"bit": "cpp",
"*.tcc": "cpp",
"chrono": "cpp",
"codecvt": "cpp",
"compare": "cpp",
"concepts": "cpp",
"exception": "cpp",
"functional": "cpp",
"iterator": "cpp",
"memory": "cpp",
"memory_resource": "cpp",
"numeric": "cpp",
"random": "cpp",
"source_location": "cpp",
"system_error": "cpp",
"type_traits": "cpp",
"utility": "cpp",
"format": "cpp",
"numbers": "cpp",
"ranges": "cpp",
"semaphore": "cpp",
"stop_token": "cpp",
"thread": "cpp",
"cfenv": "cpp",
"cinttypes": "cpp",
"typeindex": "cpp"
}
}

View File

@ -1,7 +1,13 @@
/*==================
Thie event builder both sort and build event. skip the *.to file.
===================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <cmath>
#include <stdbool.h>
#include "TFile.h"
@ -11,170 +17,220 @@
#include "TStopwatch.h"
#include "TTreeIndex.h"
#include "../mapping.h"
#include <sys/time.h> /** struct timeval, select() */
inline unsigned int getTime_us(){
unsigned int time_us;
struct timeval t1;
struct timezone tz;
gettimeofday(&t1, &tz);
time_us = (t1.tv_sec) * 1000 * 1000 + t1.tv_usec;
return time_us;
}
// #include <chrono>
// inline unsigned long long getTime_ns(){
// std::chrono::high_resolution_clock::time_point currentTime = std::chrono::high_resolution_clock::now();
// std::chrono::nanoseconds nanoseconds = std::chrono::duration_cast<std::chrono::nanoseconds>(currentTime.time_since_epoch());
// return nanoseconds.count();
// }
#include "evtReader.h"
#define MAXMULTI 100
#define BUFFERSIZE 1000000 // number of time and filePos in buffer
#define DEBUG 0
int main(int argn, char **argv){
printf("=====================================\n");
printf("=== Event Builder from *_raw.root ===\n");
printf("=== Event Builder from *.evt ===\n");
printf("=====================================\n");
if (argn != 2 && argn != 3 && argn != 4 ) {
if (argn < 6 ) {
printf("Usage :\n");
printf("%s [_raw.root File] <timeWindows> <SaveFileName>\n", argv[0]);
printf(" timeWindows : default = 100 \n");
printf(" SaveFileName : default is *.root \n");
printf("%s [timeWindows] [mapping file] [Reject Flag] [SaveFileName] [*.evt File1] [*.evt File2] ...\n", argv[0]);
printf(" timeWindows [int]: 1 unit = 10 ns \n");
printf(" mapping file path [str]: the path of mapping file. \n");
printf(" Reject Flag [int]: 0 = no rejection, 1 = reject no gamma, 2 = reject no GAGG. see mapping.h\n");
printf(" SaveFileName [str]: custom save file name \n");
return 1;
}
TString inFileName = argv[1]; // need to check name
int timeWindow = 100;
if( argn >= 3 ) timeWindow = atoi(argv[2]);
int timeWindow = atoi(argv[1]);
TString mappingFilePath = argv[2];
unsigned short rejectFlag = atoi(argv[3]);
TString outFileName = argv[4];
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;
std::vector<std::string> inFileList;
for( int i = 5; i < argn; i++) inFileList.push_back(argv[i]);
printf(" Mapping file Path : %s \n", mappingFilePath.Data());
printf(" Time window : %d ticks\n", timeWindow);
printf(" Reject Flag : %u \n", rejectFlag);
printf(" outfile Name : %s \n", outFileName.Data());
printf("--------------- Number of in files %ld \n", inFileList.size());
for( size_t i = 0; i < inFileList.size(); i++){
printf("%2ld | %s \n", i, inFileList[i].c_str());
}
TTree * tree = (TTree *) inFile->Get("tree");
printf("================================== Digesting Mapping file (to be impletmented\n");
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("================================== Creating Tree\n");
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");
UInt_t eventID = 0 ;
UInt_t multi = 0; /// this is total multipicilty for all detectors
newtree->Branch("multi", &multi, "multi/i");
newtree->Branch("evID", &eventID, "event_ID/i");
// Int_t multiCry = 0 ; /// thi is total multiplicity for all crystal
// newtree->Branch("multiCry", &multiCry, "multiplicity_crystal/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()};
UInt_t id[MAXMULTI] = {0};
Int_t e[MAXMULTI] = {-1};
ULong64_t e_t[MAXMULTI] = {0};
newtree->Branch("id", id, "id[multi]/I" );
newtree->Branch("e", e, "e[multi]/D" );
newtree->Branch("id", id, "id[multi]/i" );
newtree->Branch("e", e, "e[multi]/I" );
newtree->Branch("e_t", e_t, "e_timestamp[multi]/l");
saveFile->cd();
printf("================== Start processing....\n");
Float_t Frac = 0.05; ///Progress bar
TStopwatch StpWatch;
StpWatch.Start();
int multiOverflow = 0;
std::vector<timePos> hitList;
evtReader * reader = nullptr;
std::vector<DataBlock> event;
event.clear();
unsigned long int totalBlock;
unsigned long int blockCount = 0 ;
for( Long64_t entry = 0; entry < totnumEntry; entry++){
unsigned long long tStart = 0;
unsigned long long tEnd = 0;
/*********** 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;
unsigned int runStartTime = getTime_us();
for( size_t i = 0 ; i < inFileList.size(); i++){
reader = new evtReader(inFileList[i]);
reader->ScanNumberOfBlock();
totalBlock = reader->GetNumberOfBlock();
blockCount = 0;
do{
hitList = reader->ReadBatchPos(BUFFERSIZE, DEBUG);
blockCount += hitList.size();
printf(" %lu block %lu / %lu | %u \n", hitList.size(), blockCount, totalBlock, eventID);
if( hitList.size() == 0 ) break;
for( size_t k = 0; k < hitList.size(); k++ ){
reader->ReadBlockAtPos(hitList[k].inFilePos);
if( eventID == 0 ) tStart = reader->data->time;
if( event.size() == 0 ) {
event.push_back(*(reader->data));
if( timeWindow < 0 ) { // no event build
multi = 1;
e[0] = event[0].energy;
e_t[0] = event[0].time;
if( DEBUG ){
printf("====================== event %u, event size %u\n", eventID, multi);
printf("%6d, %12llu \n", event[0].energy, event[0].time);
}
Long64_t ev = index[entry];
saveFile->cd();
newtree->Fill();
b_ID->GetEntry(ev, 0);
b_energy->GetEntry(ev, 0);
b_energy_timestamp->GetEntry(ev, 0);
if( time0 == 0) {
time0 = energy_t;
multi = 0;
eventID ++;
event.clear();
}
timeDiff = (int) (energy_t - time0);
if( timeDiff < timeWindow ) {
continue;
if( multi > MAXMULTI ){
multiOverflow++;
}else{
id[multi] = detID;
e[multi] = energy;
e_t[multi] = energy_t;
multi ++;
if( detID < NCRYSTAL ) multiCry++;
}
if( hitList[k].time - event.front().time <= timeWindow ){
event.push_back(*(reader->data));
}else{
///---- end of event
//save event
multi = event.size();
if( DEBUG ) printf("====================== event %u, event size %u\n", eventID, multi);
for( size_t p = 0; p < multi; p++ ) {
if( DEBUG ) printf("%lu | %6d, %12llu \n", p, event[p].energy, event[p].time);
e[p] = event[p].energy;
e_t[p] = event[p].time;
}
saveFile->cd();
newtree->Fill();
eventID ++;
//clear event
event.clear();
event.push_back(*(reader->data));
}
}
}
}while(true);
tEnd = reader->data->time;
delete reader;
}
//save the last event
if( timeWindow >= 0 ){
multi = event.size();
for( size_t p = 0; p < multi; p++ ) {
if( DEBUG ) printf("%lu | %6d, %12llu \n", p, event[p].energy, event[p].time);
e[p] = event[p].energy;
e_t[p] = event[p].time;
}
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();
unsigned int runEndTime = getTime_us();
double runTime = (runEndTime - runStartTime) * 1e-6;
printf("========================================= finished.\n");
printf(" event building time = %.2f sec = %.2f min\n", runTime, runTime/60.);
printf(" total events built = %u by event builder (%llu in tree)\n", eventID , newtree->GetEntriesFast());
double tDuration_sec = (tEnd - tStart) * 1e-9;
printf(" first timestamp = %20llu ns\n", tStart);
printf(" last timestamp = %20llu ns\n", tEnd);
printf(" total data duration = %.2f sec = %.2f min\n", tDuration_sec, tDuration_sec/60.);
printf("==============> saved to %s \n", outFileName.Data());
// TMacro info;
// info.AddLine(Form("tStart= %20llu ns",tStart));
// info.AddLine(Form(" tEnd= %20llu ns",tEnd));
// info.Write("info");
saveFile->Close();
printf(" total number of event Built : %d \n", eventID);
printf(" total event has multi > %6d : %d \n", MAXMULTI, multiOverflow);
return 0;
}

View File

@ -1,68 +0,0 @@
/*==================
Thie event builder both sort and build event. skip the *.to file.
===================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cmath>
#include <stdbool.h>
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TBenchmark.h"
#include "TStopwatch.h"
#include "TTreeIndex.h"
#include "../mapping.h"
#include "evtReader.h"
#define MAXMULTI 100
#define MAXBUFFER 100 // number of hit in memory buffer
int main(int argn, char **argv){
printf("=====================================\n");
printf("=== Event Builder from *.evt ===\n");
printf("=====================================\n");
if (argn < 5 ) {
printf("Usage :\n");
printf("%s [timeWindows] [Reject Flag] [SaveFileName] [*.evt File1 [*.evt File2] ...\n", argv[0]);
printf(" timeWindows [int]: 1 unit = 10 ns \n");
printf(" Reject Flag [int]: 0 = no rejection, 1 = reject no gamma, 2 = reject no GAGG. see mapping.h\n");
printf(" SaveFileName [str]: custom save file name \n");
return 1;
}
unsigned short timeWindow = atoi(argv[1]);
unsigned short rejectFlag = atoi(argv[2]);
TString outFileName = argv[3];
std::vector<std::string> inFileList;
for( int i = 4; i < argn; i++) inFileList.push_back(argv[i]);
//========== put data into RAM buffer.
std::vector<DataBlock> hitListA, hitListB;
evtReader * reader = nullptr;
for( size_t i = 0 ; i < inFileList.size(); i++){
reader = new evtReader(inFileList[i]);
reader->ScanNumberOfBlock();
reader->ReadBatchDataAndSort(MAXBUFFER);
hitListA = reader->GetDataList();
reader->ReadBatchDataAndSort(MAXBUFFER);
hitListB = reader->GetDataList();
delete reader;
}
return 0;
}

View File

@ -1,206 +0,0 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#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] <timeWindows> <SaveFileName>\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);
}

View File

@ -1,132 +0,0 @@
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string.h>
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TMath.h"
#include "TBenchmark.h"
#include <vector>
#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);
}

File diff suppressed because it is too large Load Diff

View File

@ -17,6 +17,24 @@
#define MAX_CHANNELS_PER_BOARD 16
#define BOARD_START 2
class timePos{
public:
timePos(ULong64_t time, unsigned long int inFilePos){
this->time = time;
this->inFilePos = inFilePos;
}
ULong64_t time;
unsigned long int inFilePos;
void Print(){
printf("time: %16llu, filePos: %lu\n", time, inFilePos);
}
};
class evtReader{
public:
@ -26,6 +44,7 @@ class evtReader{
FILE * inFile;
long int inFileSize;
long int inFilePos0;
long int inFilePos;
bool endOfFile;
bool isOpened;
@ -40,8 +59,6 @@ class evtReader{
long int inFilePosPrecent[10];
Long64_t blockIDPrecent[10];
std::vector<DataBlock> dataList;
///============================================ Methods
public:
@ -65,11 +82,17 @@ class evtReader{
int ReadBlock(int opt = 0); /// 0 = default, fill data
/// 1 = no fill data
int ReadBlockAtPos(unsigned long int filePos);
void ScanNumberOfBlock();
void JumptoPrecent(int precent); ///this is offset by 1 block
void PrintStatus(int mod);
std::vector<DataBlock> ReadBatchDataAndSort(long numData);
std::vector<timePos> timePosList;
unsigned long int filePosK;
std::vector<timePos> ReadBatchPos(long batchSize, bool verbose = false);
void SortTimePos();
};
@ -82,6 +105,8 @@ evtReader::evtReader(){
inFileSize = 0;
inFilePos = 0;
inFilePos0 = 0;
filePosK = 0;
nBlock = 0;
blockID = -1;
@ -92,8 +117,7 @@ evtReader::evtReader(){
evtReader::~evtReader(){
fclose(inFile);
delete inFile;
fclose(inFile); //already delete inFile
delete data;
}
@ -104,6 +128,8 @@ evtReader::evtReader(TString inFileName){
inFileSize = 0;
inFilePos = 0;
inFilePos0 = 0;
filePosK = 0;
nBlock = 0;
blockID = -1;
@ -154,6 +180,13 @@ bool evtReader::IsEndOfFile() {
return haha > 0 ? true: false;
}
int evtReader::ReadBlockAtPos(unsigned long int filePos){
fseek( inFile, filePos, SEEK_SET);
return ReadBlock(0);
}
int evtReader::ReadBlock(int opt){
@ -162,6 +195,8 @@ int evtReader::ReadBlock(int opt){
unsigned int header[4]; ///read 4 header, unsigned int = 4 byte = 32 bits.
inFilePos0 = inFilePos;
if ( fread(header, sizeof(header), 1, inFile) != 1 ) {
endOfFile = true;
return -1;
@ -277,6 +312,8 @@ void evtReader::ScanNumberOfBlock(){
printf("scan complete: number of data Block : %ld\n", nBlock);
inFilePos = 0;
inFilePos0 = 0;
filePosK = 0;
blockID = -1;
rewind(inFile); ///back to the File begining
@ -310,23 +347,141 @@ void evtReader::PrintStatus(int mod){
}
std::vector<DataBlock> evtReader::ReadBatchDataAndSort(long numData){
std::vector<DataBlock> dataList_A;
if( endOfFile ){
return dataList_A;
void evtReader::SortTimePos(){
std::sort(timePosList.begin(), timePosList.end(), [](const timePos & a, const timePos & b) { return a.time < b.time; });
}
dataList.clear();
std::vector<timePos> evtReader::ReadBatchPos(long batchSize, bool verbose){
for( long long k = 0; k < std::min(numData, nBlock); k++){
ReadBlock();
dataList.push_back(*data);
if( verbose ) printf("&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& %s / %lu / %lu\n", __func__, filePosK, inFileSize );
fseek( inFile, filePosK, SEEK_SET);
std::vector<timePos> timePosList_A;
if( filePosK == inFileSize ){
timePosList_A = timePosList;
timePosList.clear();
return timePosList_A;
}
//sort the dataList
std::sort(dataList.begin(), dataList.end(), [](const DataBlock & a, const DataBlock & b) { return a.time < b.time; });
if( timePosList.size() == 0 ){
int res = 0;
do{
res = ReadBlock();
timePosList.push_back( timePos(data->time, inFilePos0));
}while( timePosList.size() < batchSize && res == 1);
SortTimePos();
unsigned long long t0_B = timePosList.front().time;
unsigned long long t1_B = timePosList.back().time;
if( verbose ) {
printf(" hit in memeory : %7zu | %lu | %lu \n", timePosList.size(), inFilePos, inFileSize);
printf("t0 : %15llu\n", t0_B);
printf("t1 : %15llu\n", t1_B);
}
timePosList_A = timePosList;
timePosList.clear();
}else{
timePosList_A = timePosList;
timePosList.clear();
}
if( feof(inFile) ) return timePosList_A;
int res = 0;
do{
res = ReadBlock();
timePosList.push_back(timePos(data->time, inFilePos0));
}while( timePosList.size() < batchSize && res == 1);
SortTimePos();
unsigned long long t0_B = timePosList.front().time;
unsigned long long t1_B = timePosList.back().time;
if( verbose ) {
printf(" hit in memeory : %7zu | %lu | %lu \n", timePosList.size(), inFilePos, inFileSize);
printf("t0 : %15llu\n", t0_B);
printf("t1 : %15llu\n", t1_B);
}
unsigned long long t0_A = timePosList_A.front().time;
unsigned long long t1_A = timePosList_A.back().time;
ulong ID_A = 0;
ulong ID_B = 0;
if( t0_A >= t0_B) {
printf("\033[0;31m!!!!!!!!!!!!!!!!! %s | Need to increase the batch size. \033[0m\n", __func__);
printf("present batch size : %lu \n", batchSize);
printf("timePosList_A t0_A : %15llu\n", t0_A);
printf("timePosList t0_B : %15llu\n", t0_B);
return std::vector<timePos> ();
}
if( t1_A > t0_B) { // need to sort between two hitList
if( verbose ) {
printf("############# need to sort \n");
printf("=========== sume of A + B : %zu \n", timePosList_A.size() + timePosList.size());
}
std::vector<timePos> timePosListTemp;
// find the hit that is >= t0_B, save them to timePosListTemp
for( size_t j = 0; j < timePosList_A.size() ; j++){
if( timePosList_A[j].time < t0_B ) continue;;
if( ID_A == 0 ) ID_A = j;
timePosListTemp.push_back(timePosList_A[j]);
}
// remove timePosList_A element that is >= t0_B
timePosList_A.erase(timePosList_A.begin() + ID_A, timePosList_A.end() );
// find the hit that is <= t1_A, save them to timePosListTemp
for( size_t j = 0; j < timePosList.size(); j++){
if( timePosList[j].time > t1_A ) {
break;
}
timePosListTemp.push_back(timePosList[j]);
ID_B = j + 1;
}
// remove hit elements that is <= t1_A
timePosList.erase(timePosList.begin(), timePosList.begin() + ID_B );
// sort timePosListTemp
std::sort(timePosListTemp.begin(), timePosListTemp.end(), [](const timePos& a, const timePos& b) {
return a.time < b.time;
});
if( verbose ) {
printf("----------------- ID_A : %lu, Drop\n", ID_A);
printf("----------------- ID_B : %lu, Drop\n", ID_B);
printf("=========== sume of A + B + Temp : %zu \n", timePosList_A.size() + timePosList.size() + timePosListTemp.size());
printf("----------------- refill timePosList_A \n");
}
for( size_t j = 0; j < timePosListTemp.size(); j++){
timePosList_A.push_back(timePosListTemp[j]);
}
timePosListTemp.clear();
if( verbose ) {
printf("=========== sume of A + B : %zu \n", timePosList_A.size() + timePosList.size());
printf(" A in memeory : %7zu \n", timePosList_A.size());
printf("t0 : %15llu\n", timePosList_A.front().time);
printf("t1 : %15llu\n", timePosList_A.back().time);
printf(" B in memeory : %7zu | %lu | %lu \n", timePosList.size(), inFilePos, inFileSize);
printf("t0 : %15llu\n", timePosList.front().time);
printf("t1 : %15llu\n", timePosList.back().time);
}
}
filePosK = inFilePos;
return timePosList_A;
}
#endif

View File

@ -1,11 +1,11 @@
CC=g++
#all: to2root evt2hist MergeEVT ev22txt EventBuilder pxi-time-order
all: to2root evt2hist MergeEVT ev22txt EventBuilder pxi-fsu-time-order
all: to2root evt2hist ev22txt EventBuilder pxi-fsu-time-order
#this is FSU evt to root
xia2root: ../armory/xia2root.cpp
$(CC) ../armory/xia2root.cpp -o xia2root `root-config --cflags --glibs`
# xia2root: ../armory/xia2root.cpp
# $(CC) ../armory/xia2root.cpp -o xia2root `root-config --cflags --glibs`
#xia2ev2_nopart: armory/xia2ev2_nopart.cpp
# $(CC) armory/xia2ev2_nopart.cpp -o xia2ev2_nopart
@ -15,8 +15,8 @@ to2root: ../armory/to2root.cpp ../armory/DataBlock.h ../armory/evtReader.h ../ma
$(CC) ../armory/to2root.cpp -o to2root `root-config --cflags --glibs`
#this is for online root
MergeEVT: ../armory/MergeEVT.cpp ../armory/DataBlock.h ../armory/evtReader.h ../mapping.h
$(CC) ../armory/MergeEVT.cpp -o MergeEVT `root-config --cflags --glibs`
# 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 ../armory/evtReader.h ../mapping.h

Binary file not shown.

Before

Width:  |  Height:  |  Size: 654 KiB

After

Width:  |  Height:  |  Size: 99 KiB

View File

@ -1,431 +0,0 @@
/*********************************************************/
/* 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 <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <unistd.h>
#include <time.h>
#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<MAX_ID; i++){
subevents[i] = (struct subevent *) malloc(sizeof(struct subevent)*DEF_SUB_EVENTS);
if (subevents[i] == NULL) {
printf("malloc failed\n");
return -1;
}
totmem += sizeof(struct subevent)*DEF_SUB_EVENTS;
maxevts[i] = DEF_SUB_EVENTS;
for (j=0; j<DEF_SUB_EVENTS; j++) {
subevents[i][j].data = NULL;
subevents[i][j].length = 0;
subevents[i][j].timestamp = 0;
}
}
printf("Static Memory = %ld KB (cf. MAX_ID=%d)\n", sizeof(subevents)/1024, MAX_ID);
while (1) { //main while loop
/////////
while (1) { //fill buffers until (A) maxevents or (maxevents and 2GB) is reached for any ID
//(B) EOF
//(C) auto online mode will wait for updates and break out of fill buffers for narrow time window
//read 4-byte header
if (pause == 1) {
pause = 0;
}
else {
//////////////
//auto online
while ( (fprsize - nwords*sizeof(int) < MAX_SUB_LENGTH*sizeof(int)) && online == 1) {
online = 0;
usleep(100000); //wait 0.1 seconds before checking (prevents excessive cpu usage)
//check new file size
fprpos = ftell(fpr);
fseek(fpr, 0L, SEEK_END);
fprsize = ftell(fpr);
fseek(fpr, fprpos, SEEK_SET);
//check for lock file and active PID
FPLOCK = fopen(lockfile, "r");
if (FPLOCK != NULL) {
fscanf(FPLOCK, "%d", &lockpid);
fclose(FPLOCK);
if (getpgid(lockpid) >= 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_new; j++) {
subevents[id][j].data = NULL;
subevents[id][j].length = 0;
subevents[id][j].timestamp = 0;
}
// if circular buffer is wrapped around (i.e., final is to left of intial, move data to right of initial)
if (iptr[id] + nevts[id] > evts_old) {
for (j=0; j<iptr[id] + nevts[id] - evts_old; j++) {
if (subevents[id][evts_old+j].data == NULL) {
subevents[id][evts_old+j].data = (unsigned int *) malloc(sizeof(unsigned int)*subevents[id][j].length);
if (subevents[id][evts_old+j].data == NULL) {
printf("malloc failed\n");
return -1;
}
totmem += sizeof(unsigned int)*subevents[id][j].length;
}
subevents[id][evts_old+j].length = subevents[id][j].length;
subevents[id][evts_old+j].timestamp = subevents[id][j].timestamp;
for (i=0; i<subevents[id][evts_old+j].length; i++) {
subevents[id][evts_old+j].data[i]=subevents[id][j].data[i];
}
//free data memory until it's needed again
free(subevents[id][j].data);
subevents[id][j].data = NULL;
totmem -= sizeof(unsigned int)*subevents[id][j].length;
subevents[id][j].length = 0;
subevents[id][j].timestamp = 0;
}
}
}
}
// time control of buffer filling for auto online mode (reset if initial value or large gap > 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<MAX_ID; i++){
free(subevents[i]);
totmem -= sizeof(struct subevent)*maxevts[i];
}
//print statistics last time
printf("\33[2K");
printf("Malloc (%d MB) : evts in (\x1B[34m%lld\x1B[0m) : evts out (\x1B[32m%lld\x1B[0m) : diff (\x1B[31m%lld\x1B[0m)\n", (totmem)/1024/1024, evts_tot_read, evts_tot_write, evts_tot_read-evts_tot_write);
if (outoforder > 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;
}

View File

@ -1,73 +0,0 @@
/*==================
Thie event builder both sort and build event. skip the *.to file.
===================*/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cmath>
#include <stdbool.h>
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TBenchmark.h"
#include "TStopwatch.h"
#include "TTreeIndex.h"
#include "../mapping.h"
#include "evtReader.h"
#define MAXMULTI 100
#define MAXBUFFER 100 // number of hit in memory buffer
int main(int argn, char **argv){
printf("=====================================\n");
printf("=== Event Builder from *.evt ===\n");
printf("=====================================\n");
if (argn < 5 ) {
printf("Usage :\n");
printf("%s [timeWindows] [Reject Flag] [SaveFileName] [*.evt File1 [*.evt File2] ...\n", argv[0]);
printf(" timeWindows [int]: 1 unit = 10 ns \n");
printf(" Reject Flag [int]: 0 = no rejection, 1 = reject no gamma, 2 = reject no GAGG. see mapping.h\n");
printf(" SaveFileName [str]: custom save file name \n");
return 1;
}
unsigned short timeWindow = atoi(argv[1]);
unsigned short rejectFlag = atoi(argv[2]);
TString outFileName = argv[3];
std::vector<std::string> inFileList;
for( int i = 4; i < argn; i++) inFileList.push_back(argv[i]);
//========== put data into RAM buffer.
std::vector<DataBlock> hitListA, hitListB;
evtReader * reader = nullptr;
for( size_t i = 0 ; i < inFileList.size(); i++){
reader = new evtReader(inFileList[i]);
reader->ScanNumberOfBlock();
// reader->ReadBatchDataAndSort(MAXBUFFER);
// hitListA = reader->GetDataList();
// reader->ReadBatchDataAndSort(MAXBUFFER);
// hitListB = reader->GetDataList();
delete reader;
}
return 0;
}

View File

@ -1,375 +0,0 @@
#include <stdio.h>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <string.h>
#include "TFile.h"
#include "TTree.h"
#include "TString.h"
#include "TMath.h"
#include <vector>
#define NUMDET 64 /// number of detector
#define STARTDETID 15
std::vector<std::string> SplitStr(std::string tempLine, std::string splitter, int shift = 0){
std::vector<std::string> 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<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("%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 "<<argv[0]<<" filein.evt fileout correctfilename STARTDETIDnumber maxtime start_buf# stop_buf#\n
/// Converts physics buffers (type 30) from \'file.evt\' to \'fileout.ev2\' from \'start_buf#\' until\n
/// 'stop_buf#' or whole file if no event numbers are given or \'start buf#\' = 0\n
/// STARTDETIDnum is the number of the dE ADC on a scale of 1 to n.\n
/// maxtime defaults to 100 if not specified\n
/// **Use \'NULL\' for the output file to avoid writing an ev2 file. \n
/// Note that the start and end buffer numbers are for the input where each adc \n
/// read forms another buffer.\n This version sets the time in event of the start detector to 100 and all other times relative to that.
/// If a correct file is specified, it dorrects the detectors in the file for walk.
/// It does not write out events which do not contain the start detector.\n
/// Also it only writes out events in which every ADC times in within the limits specified in the correct.tab file."<<std::endl;
exit(0);
}
printf("=======================================================\n");
printf("=== XIA evt file to CERN ROOT ===\n");
printf("=======================================================\n");
printf("The start detector number is %d\n", STARTDETID);
unsigned long long blockNum=0; /// this is total block
unsigned long long block30Num=0; /// this is number of block-30
unsigned long long eventID=0; /// this is eventID
///for( int i = 0 ; i < argn; i++) printf("............. %s \n", argv[i]);
int rawOpt = 0;
if ( argn >= 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 "<<argv[1]<<std::endl;
exit(0);
}
TString inFileName = argv[1];
TString outFileName = inFileName;
outFileName.Remove(inFileName.First('.'));
if( rawOpt == 0 ) outFileName.Append("_raw");
outFileName.Append(".root");
printf(" In file : %s \n", inFileName.Data());
printf("Out file : %s \n", outFileName.Data());
std::vector<std::vector<double>> 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;
}