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

Compare commits

merge into: rtang:master
Branches Tags
rtang:master
rtang:devel_jake
...
pull from: rtang:devel_jake
Branches Tags
rtang:devel_jake
rtang:master

1 Commits
master ... devel_jake

Author SHA1 Message Date
Jacob Davis e622f41dab Initial commit, adding all of my analysis scripts and changed how traces and cfds are handled in to2sort.cpp, datablock.cpp. 2026-06-30 15:31:35 -04:00
33 changed files with 4425 additions and 1273 deletions
Show Stats Expand all files Collapse all files

BIN
JakeStuff/DrawClovers Executable file
View File

Binary file not shown.

80
JakeStuff/DrawClovers.cpp Normal file
View File

@ -0,0 +1,80 @@
void DrawClovers(const char* filename = "outputs/Nov14_17_19.root", Bool_t drawCrystals = 1) {
TFile* file = TFile::Open(filename);
if (!file || file->IsZombie()) {
std::cerr << "Error: Could not open file " << filename << std::endl;
return;
}
const int nCanvases = 11;
const int nPerCanvas = 6;
// ROOT color codes: Blue=4, Black=1, Green=3, Red=2, Violet=6
const int colors[nPerCanvas] = {4, 1, 3, 2, 6, 30};
if (drawCrystals == true){
for (int i = 0; i < nCanvases; ++i) {
TString canvasName = Form("canvas%d", i);
TCanvas* c = new TCanvas(canvasName, canvasName, 1200, 800);
c->Divide(1, nPerCanvas);
for (int j = 0; j < nPerCanvas; ++j) {
c->cd(j + 1);
TH1F* hist = nullptr;
TString histName;
if (j < 4) {
int index = i * 4 + j;
histName = Form("CompSupCalGeEnergy%02d", index);
} else if (j==5) {
histName = Form("AddBackGeEnergy%02d", i);
}else{
histName = Form("SummedCloverCalGeEnergy%02d", i);
}
std::cout << "Trying to load histogram: " << histName << std::endl;
hist = dynamic_cast<TH1F*>(file->Get(histName));
if (hist) {
hist->SetTitle(histName); // Optional
hist->SetLineColor(colors[j]);
hist->Draw();
} else {
std::cerr << "Missing histogram: " << histName << std::endl;
}
}
c->Update();
}
}else{
const int nPerCanvas2 = 2;
const int colors2[nPerCanvas2] = {6, 30};
for (int i = 0; i < nCanvases; ++i) {
TString canvasName = Form("canvas%d", i);
TCanvas* c = new TCanvas(canvasName, canvasName, 1200, 800);
c->Divide(1, nPerCanvas2);
for (int j = 0; j < nPerCanvas2; ++j) {
c->cd(j + 1);
TH1F* hist = nullptr;
TString histName;
if (j == 0) {
histName = Form("SummedCloverCalGeEnergy%02d", i);
}else{
histName = Form("AddBackGeEnergy%02d", i);
}
std::cout << "Trying to load histogram: " << histName << std::endl;
hist = dynamic_cast<TH1F*>(file->Get(histName));
if (hist) {
hist->SetTitle(histName); // Optional
hist->SetLineColor(colors2[j]);
hist->Draw();
} else {
std::cerr << "Missing histogram: " << histName << std::endl;
}
}
c->Update();
}
}
}

6
JakeStuff/JakeToDo.txt Normal file
View File

@ -0,0 +1,6 @@
figure out binning
Make the analyzer batch runable
I need to chain a bunch of files together here I think
Or sum all of the output files at the end?
Make the entire thing batch runable
Start with .evts, run one command and end up with one root file for all run files.

BIN
JakeStuff/Jericho Executable file
View File

Binary file not shown.

1567
JakeStuff/Jericho.cpp Normal file
View File

File diff suppressed because it is too large Load Diff

19
JakeStuff/Makefile Normal file
View File

@ -0,0 +1,19 @@
# Makefile
CXX = g++
CXXFLAGS = -O2 -Wall
ROOTCFLAGS = $(shell root-config --cflags)
ROOTLIBS = $(shell root-config --libs)
TARGETS = Jericho sumOutputs
all: $(TARGETS)
Jericho: Jericho.cpp
$(CXX) $(CXXFLAGS) $(ROOTCFLAGS) -o $@ $< $(ROOTLIBS)
sumOutputs: sumOutputs.cpp
$(CXX) $(CXXFLAGS) $(ROOTCFLAGS) -o $@ $< $(ROOTLIBS)
clean:
rm -f $(TARGETS)

49
JakeStuff/batchRunJericho.sh Executable file
View File

@ -0,0 +1,49 @@
#!/bin/bash
# Usage check
if [ $# -ne 2 ]; then
echo "Usage: $0 <input_directory> <parameter_file>"
exit 1
fi
INPUT_DIR="$1"
PARAM_FILE="$2"
# Check inputs
if [ ! -d "$INPUT_DIR" ]; then
echo "Error: Directory '$INPUT_DIR' not found"
exit 1
fi
if [ ! -f "$PARAM_FILE" ]; then
echo "Error: Parameter file '$PARAM_FILE' not found"
exit 1
fi
# Loop over ROOT files
shopt -s nullglob
for rootfile in "$INPUT_DIR"/*.root; do
echo "=============================================="
echo "Processing file: $rootfile"
echo "=============================================="
# Export so runJericho.sh can see it
export inputDataFile="$rootfile"
./runJericho.sh "$PARAM_FILE"
if [ $? -ne 0 ]; then
echo "❌ Error processing $rootfile"
exit 1
fi
done
echo "✅ All files processed successfully"
source "$PARAM_FILE"
baseName1=$(basename "$inputDataFile" .root)
outPutFileName="${baseName1}_total.root"
echo "Summing outputs to $outPutFileName"
./sumOutputs $outPutDirectory $outPutFileName

62
JakeStuff/calNov2025.dat Normal file
View File

@ -0,0 +1,62 @@
0 0 0
1 0.151622084893127 0.173443485122107
2 0.149201243550541 0.231056683126326
3 0.14804661589375 0.169770867217949
4 0.154385176354815 0.27436544374234
5 0.156996807683194 0.178323431539184
6 0.151970593107449 0.149664715074664
7 0.157366291176101 0.0943410767629302
8 0.151950563195605 0.429349574565322
9 0.152799551941885 0.450677211593302
10 0.152909822993851 0.10583086021586
11 0.154785176118751 0.2176540784792
12 0.377097681276902 0.293713618342281
13 0.369625637976828 0.0450975373187248
14 0.360359104078422 -0.0549249250103685
15 0.184560138822004 0.205634627700078
16 0.129858358137908 0.288412800590095
17 0.129734152855513 0.0842495261921385
18 0.147244056087955 0.362215099497121
19 0.170507193874792 0.744583721620529
20 0 0
21 0 0
22 0 0
23 0 0
24 0.151007730564431 0.14450957772317
25 0.152332390096443 0.012267363587398
26 0.151628884271249 0.241210618081368
27 0.152223382618703 -0.120606015647013
28 0.0900641957874339 0.253536788534348
29 0.156188829257979 0.17264297616623
30 0.0889002779233048 0.427281373831534
31 0.090761895382639 0.120930296684151
32 0.150371489501676 0.191981753896812
33 0.152813923535495 0.28903652678639
34 0.151090215214167 0.160266523659516
35 0.14922170475795 0.506081115341317
36 0.155933244367767 0.214521920092466
37 0.162959149262104 0.159756894422344
38 0.168371264918135 0.151474673261532
39 0.154643360568817 0.0798072610447207
40 0 0
41 0 0
42 0 0
43 0 0
44 0.283353567841564 -0.459332495547528
45 0 0
50 0.0315226 -29.335
51 0.0315039 16.556
52 0.030819 -3.884
53 0.035185 -15.612
54 0.040979 3.1206
55 0.036493 -21.233
56 0.034558 -20.139
57 0 0
58 0.026014 70.562
59 0.031718 -162.467
60 0.027734 -1.098
61 0.026108 -0.702
62 0.026403 28.185
63 0.026531 12.483
64 0 0
65 0 0

70
JakeStuff/drawTimeDiffHists.cpp Normal file
View File

@ -0,0 +1,70 @@
// Function for GeTimeDiff histograms with 00, 01, 02, ... naming
void drawGeTimeDiffHistograms(const char* filename) {
TFile* file = TFile::Open(filename, "READ");
if (!file || file->IsZombie()) {
std::cerr << "Error: Cannot open file " << filename << std::endl;
return;
}
std::vector<TH1*> histograms;
std::vector<int> foundIndices;
// Look for histograms with GeTimeDiff00, GeTimeDiff01, etc.
for (int i = 0; i < 66; i++) {
TString histName = Form("GeTimeDiff%02d", i); // %02d gives 00, 01, 02, etc.
TH1* hist = (TH1*)file->Get(histName);
if (hist) {
histograms.push_back(hist);
foundIndices.push_back(i);
std::cout << "Found: " << histName << std::endl;
}
}
if (histograms.empty()) {
std::cerr << "No GeTimeDiff histograms found!" << std::endl;
file->Close();
return;
}
std::cout << "Found " << histograms.size() << " GeTimeDiff histograms" << std::endl;
// Create canvas with appropriate divisions
int nCols = 11;
int nRows = 6;
TCanvas* canvas = new TCanvas("canvas", "GeTimeDiff Histograms", 1920, 1080);
canvas->Divide(nCols, nRows);
// Draw histograms
for (int i = 0; i < histograms.size(); i++) {
canvas->cd(i + 1);
gPad->SetLeftMargin(0.12);
gPad->SetRightMargin(0.05);
gPad->SetTopMargin(0.1);
gPad->SetBottomMargin(0.12);
gPad->SetLogy();
histograms[i]->Draw();
// Adjust text sizes for better visibility
histograms[i]->SetTitleSize(0.06, "t");
histograms[i]->SetTitleSize(0.05, "x");
histograms[i]->SetTitleSize(0.05, "y");
histograms[i]->SetLabelSize(0.04, "x");
histograms[i]->SetLabelSize(0.04, "y");
// Add histogram index as subtitle
histograms[i]->SetTitle(Form("GeTimeDiff%02d", foundIndices[i]));
}
canvas->Update();
canvas->SaveAs("GeTimeDiff_all.pdf");
std::cout << "Canvas saved as GeTimeDiff_all.png and GeTimeDiff_all.pdf" << std::endl;
}
// Usage examples:
// root -l 'drawAllHistograms("yourfile.root")'
// root -l 'drawGeTimeDiffHistograms("yourfile.root")'

66
JakeStuff/gePromptGates.dat Normal file
View File

@ -0,0 +1,66 @@
0 -1
1 -1
2 -1
3 100
4 100
5 100
6 115
7 -1
8 -1
9 -1
10 -1
11 -1
12 125
13 100
14 125
15 125
16 125
17 125
18 150
19 150
20 150
21 150
22 150
23 150
24 125
25 125
26 200
27 200
28 200
29 200
30 -1
31 -1
32 -1
33 -1
34 -1
35 -1
36 -1
37 -1
38 -1
39 -1
40 -1
41 -1
42 100
43 175
44 130
45 150
46 140
47 150
48 130
49 140
50 130
51 160
52 100
53 90
54 150
55 150
56 150
57 150
58 150
59 150
60 -1
61 -1
62 -1
63 -1
64 -1
65 -1

66
JakeStuff/geTimingPairFile.dat Normal file
View File

@ -0,0 +1,66 @@
0 1 0
0 2 1
0 3 2
1 2 3
1 3 4
2 3 5
4 5 6
4 6 7
4 7 8
5 6 9
5 7 10
6 7 11
8 9 12
8 10 13
8 11 14
9 10 15
9 11 16
10 11 17
12 13 18
12 14 19
12 15 20
13 14 21
13 15 22
14 15 23
16 17 24
16 18 25
16 19 26
17 18 27
17 19 28
18 19 29
20 21 30
20 22 31
20 23 32
21 22 33
21 23 34
22 23 35
24 25 36
24 26 37
24 27 38
25 26 39
25 27 40
26 27 41
28 29 42
28 30 43
28 31 44
29 30 45
29 31 46
30 31 47
32 33 48
32 34 49
32 35 50
33 34 51
33 35 52
34 35 53
36 37 54
36 38 55
36 39 56
37 38 57
37 39 58
38 39 59
40 41 60
40 42 61
40 43 62
41 42 63
41 43 64
42 43 65

9
JakeStuff/jerichoTest.par Normal file
View File

@ -0,0 +1,9 @@
calFile=/home/jod23/software/XIAEventBuilder/JakeStuff/calNov2025.dat
gePromptGatesFile=/home/jod23/software/XIAEventBuilder/JakeStuff/gePromptGates.dat
outPutDirectory=/home/jod23/
trinABCoin=0
addbackDiagonals=0
alphaCutsFileName=/home/jod23/software/XIAEventBuilder/JakeStuff/alphaCuts.root
geTimingPairFile=/home/jod23/software/XIAEventBuilder/JakeStuff/geTimingPairFile.dat
bgoTimeGateNS=250
labrTimingPairFile=/home/jod23/software/XIAEventBuilder/JakeStuff/labrTimingPairFile.dat

120
JakeStuff/labrTimingPairFile.dat Normal file
View File

@ -0,0 +1,120 @@
0 1 0
0 2 1
0 3 2
0 4 3
0 5 4
0 6 5
0 7 6
0 8 7
0 9 8
0 10 9
0 11 10
0 12 11
0 13 12
0 14 13
0 15 14
1 2 15
1 3 16
1 4 17
1 5 18
1 6 19
1 7 20
1 8 21
1 9 22
1 10 23
1 11 24
1 12 25
1 13 26
1 14 27
1 15 28
2 3 29
2 4 30
2 5 31
2 6 32
2 7 33
2 8 34
2 9 35
2 10 36
2 11 37
2 12 38
2 13 39
2 14 40
2 15 41
3 4 42
3 5 43
3 6 44
3 7 45
3 8 46
3 9 47
3 10 48
3 11 49
3 12 50
3 13 51
3 14 52
3 15 53
4 5 54
4 6 55
4 7 56
4 8 57
4 9 58
4 10 59
4 11 60
4 12 61
4 13 62
4 14 63
4 15 64
5 6 65
5 7 66
5 8 67
5 9 68
5 10 69
5 11 70
5 12 71
5 13 72
5 14 73
5 15 74
6 7 75
6 8 76
6 9 77
6 10 78
6 11 79
6 12 80
6 13 81
6 14 82
6 15 83
7 8 84
7 9 85
7 10 86
7 11 87
7 12 88
7 13 89
7 14 90
7 15 91
8 9 92
8 10 93
8 11 94
8 12 95
8 13 96
8 14 97
8 15 98
9 10 99
9 11 100
9 12 101
9 13 102
9 14 103
9 15 104
10 11 105
10 12 106
10 13 107
10 14 108
10 15 109
11 12 110
11 13 111
11 14 112
11 15 113
12 13 114
12 14 115
12 15 116
13 14 117
13 15 118
14 15 119

53
JakeStuff/mappingNov2025.h Normal file
View File

@ -0,0 +1,53 @@
/************************************
Clover : 0 - 99
BGO : 100 - 199
GAGG A : 200 - 299
GAGG B : 300 - 399
ZERO DEGREE : 400 - 499
BEGE : 44 - 45
LaBr3 : 500 - 599
* *********************************/
#ifndef MAPPING
#define MAPPING
//==================== mapping
#define NCLOVER 11
#define NCRYSTAL NCLOVER*4
#define NBGO NCLOVER
#define NGAGG 26
#define NZEROGAGG 2 ///NZERO is used
//Help for people looking at this mapping:
//rows are pixie slots, colomns correspond to channels. Key for the first row:
//clover1_blue, clover1_black, clover1_green, clover1_red BGO1, clover2_blue, clover2_black, clover2_green, clover2_red, BGO2, clover3_blue, clover3_black, clover3_green, clover3_red, BGO3, EMPTY_CHANNEL
//This also does not take the missing detectors into account, it assumes a detector in each channel.
int mapping[208] ={
//***************** <-- load indicator for EventBuidler
//-0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0, 1, 2, 3, 100, 4, 5, 6, 7, 101, 8, 9, 10, 11, 102, 44, //slot 2
12, 13, 14, 15, 103, 16, 17, 18, 19, 104, 20, 21, 22, 23, 105, 45, //slot 3
24, 25, 26, 27, 106, 28, 29, 30, 31, 107, 32, 33, 34, 35, 108, -1, //slot 4
36, 37, 38, 39, 109, 40, 41, 42, 43, 110, 400, 401, -1, -1, -1, -1, //slot 5
200, 300, 201, 301, 202, 302, 203, 303, 204, 304, 205, 305, 206, 306, 207, 307, //slot 6 Trinity ring 1 (A and B channels alternating) crystal 1 - 8
208, 308, 209, 309, 210, 310, 211, 311, 212, 312, 213, 313, 214, 314, 215, 315, //slot 7 Ring 2, A and B alternating crystal 1 - 8
216, 316, 217, 317, 218, 318, 219, 319, 220, 320, 221, 321, 222, 322, 223, 323, //slot 8 Ring 2 first 4 channels ( crystals 9 - 10 A and B ) the rest is Ring 3 crystals 1 - 6 alternating A and B
224, 324, 225, 325, 226, 326, 227, 327, 228, 328, 229, 329, 230, 330, 231, 331, //slot 9 Ring 3 crystal 7 - 14 A and B alternating
232, 332, 233, 333, 234, 334, 235, 335, 236, 336, 237, 337, 238, 338, 239, 339, //slot 10 Ring 4 crystals 1 - 8
240, 340, 241, 341, 242, 342, 243, 343, 244, 344, 245, 345, 246, 346, 247, 347, //slot 11 Ring 4 crystals 9 - 16
248, 348, 249, 349, 250, 350, 251, 351, 252, 352, 253, 353, 254, 354, 255, 355, //slot 12 Ring 5 crystals 1 - 8
256, 356, 257, 357, 258, 358, 259, 359, 260, 360, 261, 361, 262, 362, 263, 363, //slot 13 Ring 5 crystals 9 -16
500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, //slot 14 LaBr3
//**************** <-- end of mapping indicator EventBuidler
};
//200- 209 GAGG 2A
//210- 225 GAGG 4A
//250- 259 GAGG 2B
//260- 275 GAGG 4B
#endif

31
JakeStuff/runJericho.sh Executable file
View File

@ -0,0 +1,31 @@
#!/bin/bash
if [ $# -ne 1 ]; then
echo "Usage: $0 <parameter file>"
exit 1
fi
PARAM_FILE="$1"
if [ ! -f "$PARAM_FILE" ]; then
echo "Error: File '$PARAM_FILE' not found!"
exit 1
fi
source "$PARAM_FILE"
# Construct output filename from input
baseName=$(basename "$inputDataFile" .root)
outPutFileName="${baseName}_output.root"
./Jericho "$inputDataFile" \
"$calFile" \
"$gePromptGatesFile" \
"$outPutDirectory" \
"$outPutFileName" \
"$trinABCoin" \
"$addbackDiagonals" \
"$alphaCutsFileName" \
"$geTimingPairFile" \
"$bgoTimeGateNS" \
"$labrTimingPairFile"

16
JakeStuff/saveCuts.cpp Normal file
View File

@ -0,0 +1,16 @@
void saveAllCuts(const char* filename = "alphaCuts.root") {
TFile* f = new TFile(filename, "UPDATE");
TIter next(gROOT->GetListOfSpecials());
TObject* obj;
while ((obj = next())) {
if (obj->InheritsFrom("TCutG")) {
obj->Write("", TObject::kOverwrite);
std::cout << "Saved: " << obj->GetName() << std::endl;
}
}
f->Close();
// Return to the originally attached file
if (_file0) _file0->cd();
}

BIN
JakeStuff/sumOutputs Executable file
View File

Binary file not shown.

116
JakeStuff/sumOutputs.cpp Executable file
View File

@ -0,0 +1,116 @@
#include <TFile.h>
#include <TKey.h>
#include <TDirectory.h>
#include <TSystem.h>
#include <TH1.h>
#include <iostream>
#include <vector>
#include <string>
#include <iomanip> // for setw
void AddDir(TDirectory* src, TDirectory* dest)
{
TIter next(src->GetListOfKeys());
TKey* key;
while ((key = (TKey*)next())) {
TObject* obj = key->ReadObj();
if (obj->InheritsFrom("TDirectory")) {
TDirectory* srcSub = (TDirectory*)obj;
TDirectory* destSub =
dest->GetDirectory(srcSub->GetName());
if (!destSub)
destSub = dest->mkdir(srcSub->GetName());
AddDir(srcSub, destSub);
}
else if (obj->InheritsFrom("TH1")) {
TH1* h = (TH1*)obj;
TH1* hsum =
(TH1*)dest->Get(h->GetName());
if (!hsum) {
hsum = (TH1*)h->Clone();
hsum->SetDirectory(dest);
hsum->Sumw2(kTRUE);
} else {
hsum->Add(h);
}
hsum->SetOption("HIST");
}
}
}
// Simple text-based progress bar
void printProgress(int current, int total, int barWidth = 50)
{
float progress = float(current) / total;
int pos = int(barWidth * progress);
std::cout << "[";
for (int i = 0; i < barWidth; ++i) {
if (i < pos) std::cout << "=";
else if (i == pos) std::cout << ">";
else std::cout << " ";
}
std::cout << "] " << int(progress * 100.0) << "%\r";
std::cout.flush();
}
int main(int argc, char** argv)
{
if (argc != 3) {
std::cerr << "Usage: " << argv[0]
<< " <directory> <output_file>\n";
return 1;
}
TString dir = argv[1];
TString outName = argv[2]; // <-- now from command-line
void* dirp = gSystem->OpenDirectory(dir);
if (!dirp) {
std::cerr << "Cannot open directory "
<< dir << std::endl;
return 1;
}
// Gather all ROOT files in a vector
std::vector<TString> files;
const char* entry;
while ((entry = gSystem->GetDirEntry(dirp))) {
TString fname(entry);
if (!fname.EndsWith(".root")) continue;
if (fname == outName) continue; // avoid overwriting output
files.push_back(fname);
}
gSystem->FreeDirectory(dirp);
if (files.empty()) {
std::cerr << "No ROOT files found in directory." << std::endl;
return 1;
}
TFile* fout = new TFile(dir + "/" + outName, "RECREATE");
// Loop over files with progress bar
int totalFiles = files.size();
for (int i = 0; i < totalFiles; ++i) {
TString full = dir + "/" + files[i];
TFile f(full, "READ");
if (!f.IsZombie())
AddDir(&f, fout);
f.Close();
printProgress(i + 1, totalFiles);
}
fout->Write();
fout->Close();
std::cout << std::endl << "Output written to " << dir + "/" + outName << std::endl;
return 0;
}

9
README.md
View File

@ -9,6 +9,7 @@ Codes need to be compiled are in armory/
# Function of Programs # Function of Programs
![Alt text](armory/programs_explain.png?raw=true "program explained") ![Alt text](armory/programs_explain.png?raw=true "program explained")
# armory/DataBlock.h # armory/DataBlock.h
this is the source file for the class DataBlock, it stored all information from a single data block from pixie16 output. this is the source file for the class DataBlock, it stored all information from a single data block from pixie16 output.
@ -17,8 +18,11 @@ this is the source file for the class evtReader.
It read the *.evt file (which is same as pixie16 output) and convert each measurement (or data block) from byte into meaningful data and use DataBlock class to store the information. It read the *.evt file (which is same as pixie16 output) and convert each measurement (or data block) from byte into meaningful data and use DataBlock class to store the information.
It can also scan the evt file. It can also scan the evt file.
# armory/MergeEVT
this merges all evt files into *_raw.root.
# armory/EventBuilder # armory/EventBuilder
This builds events from *.evt files to *.root file. It by pass the *.to file. the default BUFFERSIZE is 1 million datablock. this builds events from *_raw,root to *.root file
# armory/evt2hist # armory/evt2hist
this processes evt file to hstograms. this processes evt file to hstograms.
@ -30,6 +34,9 @@ this build events from *.evt.to files to *.root file (need to check the compacta
this sorting the time from *evt file to *.evt.to.fsu.XXX, where XXX is the time window. this sorting the time from *evt file to *.evt.to.fsu.XXX, where XXX is the time window.
It will search data within XXX time window, if non of the data is from clover, discard. It will search data within XXX time window, if non of the data is from clover, discard.
# armory/xia2root
this is old evt to root for custom pixie DAQ.
# Analyzer.C/h # Analyzer.C/h
this is a TSelector for analysis the *.root file this is a TSelector for analysis the *.root file

8
armory/AutoFit.C
View File

@ -15,10 +15,6 @@
#include <TColor.h> #include <TColor.h>
#include <TSpectrum.h> #include <TSpectrum.h>
#include <TMath.h> #include <TMath.h>
#include <TRandom.h>
#include <TLatex.h>
#include <TH1.h>
#include <TMarker.h>
#include <vector> #include <vector>
void showFitMethod(){ void showFitMethod(){
@ -210,7 +206,7 @@ void GoodnessofFit(TH1F * hist, TF1 * fit){
double x = hist->GetBinCenter(i); double x = hist->GetBinCenter(i);
double ybar = fit->Eval(x); double ybar = fit->Eval(x);
ysq += y*y; ysq += y*y;
mean += y; mean + y;
SSR += (y - ybar)*(y-ybar); SSR += (y - ybar)*(y-ybar);
chisq += (y - ybar)*(y-ybar)/e/e; chisq += (y - ybar)*(y-ybar)/e/e;
@ -1955,7 +1951,7 @@ void saveFitPara(TString fileName = "AutoFit_para.txt"){
fprintf(file_out, "# for n-Gauss fit, can use \"#\", or \"//\" to comment out whole line\n"); fprintf(file_out, "# for n-Gauss fit, can use \"#\", or \"//\" to comment out whole line\n");
fprintf(file_out, "# peak low high fixed? sigma_Max fixed? hight\n"); fprintf(file_out, "# peak low high fixed? sigma_Max fixed? hight\n");
for( int i = 0 ; i < (int) xPeakList.size() ; i++){ for( int i = 0 ; i < xPeakList.size() ; i++){
fprintf(file_out, "%.3f %.3f %.3f 0 %.3f 0 %.0f\n", fprintf(file_out, "%.3f %.3f %.3f 0 %.3f 0 %.0f\n",
xPeakList[i], xPeakList[i],
xPeakList[i] - 5*sigma[i], xPeakList[i] - 5*sigma[i],

4
armory/DataBlock.h
View File

@ -77,7 +77,7 @@ public:
} }
void Print(bool printTrace = true){ void Print(){
printf("============== eventID : %llu\n", eventID); printf("============== eventID : %llu\n", eventID);
printf("Crate: %d, Slot: %d, Ch: %d \n", crate, slot, ch); printf("Crate: %d, Slot: %d, Ch: %d \n", crate, slot, ch);
printf("HeaderLength: %d, Event Length: %d, energy: %d, timeStamp: %llu\n", headerLength, eventLength, energy, time); printf("HeaderLength: %d, Event Length: %d, energy: %d, timeStamp: %llu\n", headerLength, eventLength, energy, time);
@ -92,7 +92,7 @@ public:
printf(" QDCsum : \n"); printf(" QDCsum : \n");
for( int i = 0; i < 8; i++) printf(" %-10d\n", QDCsum[i]); for( int i = 0; i < 8; i++) printf(" %-10d\n", QDCsum[i]);
} }
if( printTrace && eventLength > headerLength ){ if( eventLength > headerLength ){
printf(" trace:\n"); printf(" trace:\n");
for( int i = 0 ; i < trace_length ; i++)printf("%3d| %-10d\n",i, trace[i]); for( int i = 0 ; i < trace_length ; i++)printf("%3d| %-10d\n",i, trace[i]);
} }

438
armory/EventBuilder.cpp
View File

@ -1,21 +1,8 @@
/*==================
Thie event builder both sort and build event. skip the *.to file.
===================*/
#include <stdio.h> #include <stdio.h>
#include <stdlib.h> #include <stdlib.h>
#include <string.h> #include <string.h>
#include <cmath> #include <math.h>
#include <stdbool.h> #include <stdbool.h>
#include <iostream>
#include <fstream>
#include <sys/time.h> /** struct timeval, select() */
#include <chrono>
#include <sstream>
#include <algorithm> // for std::remove_if
#include <cctype> // for std::isspace
#include "TFile.h" #include "TFile.h"
#include "TTree.h" #include "TTree.h"
@ -24,353 +11,170 @@ Thie event builder both sort and build event. skip the *.to file.
#include "TStopwatch.h" #include "TStopwatch.h"
#include "TTreeIndex.h" #include "TTreeIndex.h"
std::string trimSpaces(const std::string& str); #include "../mapping.h"
std::vector<std::string> split(const std::string& str, char delimiter);
unsigned int getTime_us();
unsigned long long getTime_ns();
#include "evtReader.h"
#define MAXMULTI 100 #define MAXMULTI 100
#define BUFFERSIZE 1000000 // number of time and filePos in buffer
#define DEBUG 0
int main(int argn, char **argv){ int main(int argn, char **argv){
printf("=====================================\n"); printf("=====================================\n");
printf("=== Event Builder from *.evt ===\n"); printf("=== Event Builder from *_raw.root ===\n");
printf("=====================================\n"); printf("=====================================\n");
if (argn < 6 ) { if (argn != 2 && argn != 3 && argn != 4 ) {
printf("Usage :\n"); printf("Usage :\n");
printf("%s [timeWindows] [mapping file] [Reject Flag] [QDC Flag] [SaveFileName] [*.evt File1] [*.evt File2] ...\n", argv[0]); printf("%s [_raw.root File] <timeWindows> <SaveFileName>\n", argv[0]);
printf(" timeWindows [int]: 1 unit = 10 ns \n"); printf(" timeWindows : default = 100 \n");
printf(" mapping file path [str]: the path of mapping file. \n"); printf(" SaveFileName : default is *.root \n");
printf(" Reject Flag [int]: 0 = no rejection. see mapping.h\n");
printf(" 1 = reject BGO\n");
printf(" 2 = reject no gamma\n");
printf(" 4 = reject no GAGG\n");
printf(" 8 = reject zero energy data-point\n");
printf(" 3 = reject BGO + no gamma, 5 = reject BGO + no GAGG, etc.\n");
printf(" QDC Flag [int]: 0 = no qdc, 1 = with qdc\n");
printf(" SaveFileName [str]: custom save file name \n");
return 1; return 1;
} }
int timeWindow = atoi(argv[1]); TString inFileName = argv[1]; // need to check name
TString mappingFilePath = argv[2]; int timeWindow = 100;
unsigned short rejectFlag = atoi(argv[3]); if( argn >= 3 ) timeWindow = atoi(argv[2]);
unsigned short qdcFlag = atoi(argv[4]);
TString outFileName = argv[5];
std::vector<std::string> inFileList; printf(">>> Opening input %s \n", inFileName.Data());
for( int i = 6; i < argn; i++) inFileList.push_back(argv[i]); TFile * inFile = new TFile(inFileName, "READ");
if( inFile->IsOpen() == false ) {
printf(" Mapping file Path : %s \n", mappingFilePath.Data()); printf("!!!! cannot open file %s \n", inFileName.Data());
printf(" Time window : %d ticks\n", timeWindow); return 0;
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());
} }
if( rejectFlag > 0 ) printf("================================== Rejection Filter Conditions\n"); TTree * tree = (TTree *) inFile->Get("tree");
if( ( rejectFlag & 0x1 ) ) printf("\033[31m !!!! Reject event with BGO !!!! \033[0m\n");
if( ( rejectFlag & 0x2 ) ) printf("\033[31m !!!! Reject event w/o Clover !!!! \033[0m\n");
if( ( rejectFlag & 0x4 ) ) printf("\033[31m !!!! Reject event w/o GAGG !!!! \033[0m\n");
if( ( rejectFlag & 0x8 ) ) printf("\033[31m !!!! Reject Zero-Energy Data Point !!!! \033[0m\n");
if( rejectFlag > 0 ) printf(" Rejection filter does not apply to the last event when timeWindow >= 0 \n");
printf("================================== Digesting Mapping file\n");
std::ifstream mapFile( mappingFilePath.Data() ); Long64_t evID;
UShort_t detID;
UShort_t energy;
ULong64_t energy_t;
int mapping[16*12]; // fixed to be 12 digitizers TBranch *b_data_ID; //!
for( int i = 0; i < 16*12 ; i++ ) mapping[i] = -1; TBranch *b_ID; //!
TBranch *b_energy; //!
TBranch *b_energy_timestamp; //!
if( !mapFile.is_open() ){ tree->SetBranchAddress("evID", &evID, &b_data_ID);
printf("Cannot open mapping file : %s. Skip. \n", mappingFilePath.Data()); tree->SetBranchAddress("id", &detID, &b_ID);
}else{ tree->SetBranchAddress("e", &energy, &b_energy);
tree->SetBranchAddress("e_t", &energy_t, &b_energy_timestamp);
int index = 0; Long64_t totnumEntry = tree->GetEntries();
bool startMap = false;
std::string line;
while( std::getline(mapFile, line) ){
// printf("|%s|\n", line.c_str());
if( line.find("//-") != std::string::npos ) continue;
if( line.find("<--") != std::string::npos ){
startMap = true;
continue;
}
if( startMap && line.find("<--") != std::string::npos ) break; printf(" total Entry : %lld \n", totnumEntry);
if( startMap ){
std::vector<std::string> list = split(line, ',');
for( size_t k = 0; k < list.size(); k ++ ){
if( list[k].find("//") != std::string::npos || k >= 16 ) continue;
// printf("%ld | %s \n", k, list[k].c_str());
mapping[index] = atoi(list[k].c_str());
index ++;
}
} printf(" event Build window: %d tick = %d nsec \n", timeWindow, timeWindow * 10);
}
mapFile.close();
}
//---- print mapping; printf(">>> Buidling Index using the timestamp\n");
for( int i = 0; i < 16*11; i++){ tree->BuildIndex("e_t");
if( i % 16 == 0 ) printf("Mod-%02d | ", i/16); TTreeIndex *in = (TTreeIndex*) tree->GetTreeIndex();
if( mapping[i] < 0 ) printf("%4d,", mapping[i]); Long64_t * index = in->GetIndex();
if( 0 <= mapping[i] && mapping[i] < 100 ) printf("\033[31m%4d\033[0m,", mapping[i]);
if( 100 <= mapping[i] && mapping[i] < 200 ) printf("\033[34m%4d\033[0m,", mapping[i]);
if( 200 <= mapping[i] && mapping[i] < 300 ) printf("\033[32m%4d\033[0m,", mapping[i]);
if( 300 <= mapping[i] && mapping[i] < 400 ) printf("\033[35m%4d\033[0m,", mapping[i]);
if( i % 16 == 15 ) printf("\n");
}
ULong64_t time0; //time-0 for each event
int timeDiff;
printf("================================== Creating Tree\n"); 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(">>> Create output tree\n"); printf(">>> Create output tree\n");
TFile * saveFile = new TFile(outFileName, "recreate"); TFile * saveFile = new TFile(outFileName, "recreate");
saveFile->cd();
TTree * newtree = new TTree("tree", outFileName); TTree * newtree = new TTree("tree", outFileName);
UInt_t eventID = 0 ; Int_t eventID = 0 ;
UInt_t multi = 0; /// this is total multipicilty for all detectors Int_t multi = 0; /// this is total multipicilty for all detectors
newtree->Branch("multi", &multi, "multi/i"); 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");
UInt_t id[MAXMULTI] = {0}; newtree->Branch("evID", &eventID, "event_ID/l");
Int_t e[MAXMULTI] = {-1};
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()};
ULong64_t e_t[MAXMULTI] = {0}; ULong64_t e_t[MAXMULTI] = {0};
UInt_t qdc[MAXMULTI][8] = {0}; newtree->Branch("id", id, "id[multi]/I" );
newtree->Branch("id", id, "id[multi]/i" ); newtree->Branch("e", e, "e[multi]/D" );
newtree->Branch("e", e, "e[multi]/I" );
newtree->Branch("e_t", e_t, "e_timestamp[multi]/l"); newtree->Branch("e_t", e_t, "e_timestamp[multi]/l");
if( qdcFlag ) newtree->Branch("qdc", qdc, "qdc[multi][8]/I");
saveFile->cd();
printf("================== Start processing....\n"); printf("================== Start processing....\n");
Float_t Frac = 0.05; ///Progress bar Float_t Frac = 0.05; ///Progress bar
TStopwatch StpWatch; TStopwatch StpWatch;
StpWatch.Start(); StpWatch.Start();
std::vector<timePos> hitList; int multiOverflow = 0;
evtReader * reader = nullptr;
std::vector<DataBlock> event;
event.clear();
unsigned long int totalBlock;
unsigned long int blockCount = 0 ;
unsigned long long tStart = 0; for( Long64_t entry = 0; entry < totnumEntry; entry++){
unsigned long long tEnd = 0;
unsigned int runStartTime = getTime_us(); /*********** Progress Bar ******************************************/
if (entry>totnumEntry*Frac-1) {
for( size_t i = 0 ; i < inFileList.size(); i++){ TString msg; msg.Form("%llu", totnumEntry/1000);
int len = msg.Sizeof();
reader = new evtReader(inFileList[i]); printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n",
reader->ScanNumberOfBlock(); Frac*100, len, entry/1000,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac);
totalBlock = reader->GetNumberOfBlock(); StpWatch.Start(kFALSE);
blockCount = 0; Frac+=0.1;
}
do{
Long64_t ev = index[entry];
hitList = reader->ReadBatchPos(BUFFERSIZE, DEBUG);
blockCount += hitList.size(); b_ID->GetEntry(ev, 0);
// std::cout << "please wait....."; b_energy->GetEntry(ev, 0);
printf("File-%ld %10lu block %10lu / %lu [%4.1f%%] | %u \r", i, hitList.size(), blockCount, totalBlock, blockCount*100./totalBlock, eventID); b_energy_timestamp->GetEntry(ev, 0);
fflush(stdout);
if( hitList.size() == 0 ) break; if( time0 == 0) {
time0 = energy_t;
for( size_t k = 0; k < hitList.size(); k++ ){ multi = 0;
}
reader->ReadBlockAtPos(hitList[k].inFilePos); timeDiff = (int) (energy_t - time0);
if( eventID == 0 ) tStart = reader->data->time; if( timeDiff < timeWindow ) {
if( event.size() == 0 ) { if( multi > MAXMULTI ){
multiOverflow++;
event.push_back(*(reader->data)); }else{
id[multi] = detID;
if( timeWindow < 0 ) { // no event build e[multi] = energy;
multi = 1; e_t[multi] = energy_t;
int index = event[0].crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (event[0].slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + event[0].ch; multi ++;
id[0] = mapping[index]; if( detID < NCRYSTAL ) multiCry++;
e[0] = event[0].energy; }
e_t[0] = event[0].time; }else{
///---- end of event
if( qdcFlag ) for( int i = 0; i < 8; i++) qdc[0][i] = event[0].QDCsum[i]; saveFile->cd();
newtree->Fill();
if( DEBUG ){ eventID ++;
printf("====================== event %u, event size %u\n", eventID, multi);
printf("%6d, %12llu \n", event[0].energy, event[0].time); ///---- clear data
} for( int i = 0; i < MAXMULTI ; i ++){
id[i] = 0;
if( (rejectFlag & 0x8 ) && e[0] == 0 ) { e[i] = TMath::QuietNaN();
event.clear(); e_t[i] = 0;
continue; }
} multi = 0;
multiCry = 0;
saveFile->cd();
newtree->Fill();
eventID ++; /// fill 1st data of an event
time0 = energy_t;
event.clear(); timeDiff = 0;
}
id[multi] = detID;
continue; e[multi] = energy;
e_t[multi] = energy_t;
}else{ multi++;
if( detID < NCRYSTAL ) multiCry++;
if( hitList[k].time - event.front().time <= timeWindow ){
event.push_back(*(reader->data));
}else{
//save event
if( DEBUG ) printf("====================== event %u, event size %lu\n", eventID, event.size());
int nBGO = 0;
int nClover = 0;
int nGagg = 0;
int count = 0 ;
for( size_t p = 0; p < event.size(); p++ ) {
if( (rejectFlag & 0x8 ) && event[p].energy == 0 ) continue;
e[count] = event[p].energy;
e_t[count] = event[p].time;
int index = event[p].crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (event[p].slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + event[p].ch;
id[count] = mapping[index];
if( qdcFlag ) for( int i = 0; i < 8; i++) qdc[count][i] = event[p].QDCsum[i];
if( DEBUG ) printf("%u | %3d, %6d, %12llu \n", count, id[count], e[count], e_t[count]);
if( 0 <= id[count] && id[count] < 100 ) nClover ++;
if( 100 <= id[count] && id[count] < 200 ) nBGO ++;
if( 200 <= id[count] && id[count] < 400 ) nGagg ++;
count ++;
}
multi = count;
if( DEBUG ) printf("nBGO %d, nClover %d, nGagg %d \n", nBGO, nClover, nGagg);
if( (rejectFlag & 0x1) && nBGO > 0 ) {
event.clear();
event.push_back(*(reader->data));
continue;
}
if( (rejectFlag & 0x2) && nClover == 0 ) {
event.clear();
event.push_back(*(reader->data));
continue;
}
if( (rejectFlag & 0x4) && nGagg == 0 ) {
event.clear();
event.push_back(*(reader->data));
continue;
}
if( multi > 0 ){
if( DEBUG ) printf("---------------> fill tree\n");
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);
int index = event[p].crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (event[p].slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + event[p].ch;
id[p] = mapping[index];
e[p] = event[p].energy;
e_t[p] = event[p].time;
if( qdcFlag ) for( int i = 0; i < 8; i++) qdc[p][i] = event[p].QDCsum[i];
} }
saveFile->cd();
newtree->Fill();
eventID ++;
} }
printf("============================== finished.\n");
saveFile->cd();
newtree->Write(); 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(); saveFile->Close();
printf(" total number of event Built : %d \n", eventID);
return 0; printf(" total event has multi > %6d : %d \n", MAXMULTI, multiOverflow);
}
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;
}
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();
}
std::string trimSpaces(const std::string& str) {
std::string trimmed = str;
trimmed.erase(std::remove_if(trimmed.begin(), trimmed.end(), ::isspace), trimmed.end());
return trimmed;
}
std::vector<std::string> split(const std::string& str, char delimiter) {
std::vector<std::string> tokens;
std::stringstream ss(str);
std::string token;
while (std::getline(ss, token, delimiter)) {
tokens.push_back(trimSpaces(token));
}
return tokens;
} }

206
armory/EventBuilder_seperated.cpp Normal file
View File

@ -0,0 +1,206 @@
#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);
}

132
armory/MergeEVT.cpp Normal file
View File

@ -0,0 +1,132 @@
#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);
}

1468
armory/XIAEventBuilder.excalidraw
View File

File diff suppressed because it is too large Load Diff

182
armory/evtReader.h
View File

@ -17,24 +17,6 @@
#define MAX_CHANNELS_PER_BOARD 16 #define MAX_CHANNELS_PER_BOARD 16
#define BOARD_START 2 #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{ class evtReader{
public: public:
@ -44,7 +26,6 @@ class evtReader{
FILE * inFile; FILE * inFile;
long int inFileSize; long int inFileSize;
long int inFilePos0;
long int inFilePos; long int inFilePos;
bool endOfFile; bool endOfFile;
bool isOpened; bool isOpened;
@ -82,18 +63,10 @@ class evtReader{
int ReadBlock(int opt = 0); /// 0 = default, fill data int ReadBlock(int opt = 0); /// 0 = default, fill data
/// 1 = no fill data /// 1 = no fill data
int ReadBlockAtPos(unsigned long int filePos);
void ScanNumberOfBlock(); void ScanNumberOfBlock();
void JumptoPrecent(int precent); ///this is offset by 1 block void JumptoPrecent(int precent); ///this is offset by 1 block
void PrintStatus(int mod); void PrintStatus(int mod);
std::vector<timePos> timePosList;
unsigned long int filePosK;
std::vector<timePos> ReadBatchPos(long batchSize, bool verbose = false);
void SortTimePos();
}; };
@ -105,8 +78,6 @@ evtReader::evtReader(){
inFileSize = 0; inFileSize = 0;
inFilePos = 0; inFilePos = 0;
inFilePos0 = 0;
filePosK = 0;
nBlock = 0; nBlock = 0;
blockID = -1; blockID = -1;
@ -117,7 +88,8 @@ evtReader::evtReader(){
evtReader::~evtReader(){ evtReader::~evtReader(){
fclose(inFile); //already delete inFile fclose(inFile);
delete inFile;
delete data; delete data;
} }
@ -128,8 +100,6 @@ evtReader::evtReader(TString inFileName){
inFileSize = 0; inFileSize = 0;
inFilePos = 0; inFilePos = 0;
inFilePos0 = 0;
filePosK = 0;
nBlock = 0; nBlock = 0;
blockID = -1; blockID = -1;
@ -180,13 +150,6 @@ bool evtReader::IsEndOfFile() {
return haha > 0 ? true: false; 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){ int evtReader::ReadBlock(int opt){
@ -195,8 +158,6 @@ int evtReader::ReadBlock(int opt){
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.
inFilePos0 = inFilePos;
if ( fread(header, sizeof(header), 1, inFile) != 1 ) { if ( fread(header, sizeof(header), 1, inFile) != 1 ) {
endOfFile = true; endOfFile = true;
return -1; return -1;
@ -312,8 +273,6 @@ void evtReader::ScanNumberOfBlock(){
printf("scan complete: number of data Block : %ld\n", nBlock); printf("scan complete: number of data Block : %ld\n", nBlock);
inFilePos = 0; inFilePos = 0;
inFilePos0 = 0;
filePosK = 0;
blockID = -1; blockID = -1;
rewind(inFile); ///back to the File begining rewind(inFile); ///back to the File begining
@ -347,141 +306,4 @@ void evtReader::PrintStatus(int mod){
} }
void evtReader::SortTimePos(){
std::sort(timePosList.begin(), timePosList.end(), [](const timePos & a, const timePos & b) { return a.time < b.time; });
}
std::vector<timePos> evtReader::ReadBatchPos(long batchSize, bool verbose){
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;
}
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 #endif

23
armory/makefile
View File

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

BIN
armory/programs_explain.png
View File

Binary file not shown.
Side by Side Swipe Overlay

Before

Width:  |  Height:  |  Size: 99 KiB

After

Width:  |  Height:  |  Size: 654 KiB

4
armory/pxi-fsu-time-order.cpp
View File

@ -449,11 +449,11 @@ int main(int argc, char **argv) {
}else{ }else{
if( FillFlag ){ if( FillFlag ){
if ( count < debugCount) printf("----------------- filled %d\n", count); if ( count < debugCount) printf("----------------- filled \n");
for( i = 0; i < nFill; i++) fwrite(fillevents[i].data, sizeof(unsigned int)*fillevents[i].length, 1, fpw); for( i = 0; i < nFill; i++) fwrite(fillevents[i].data, sizeof(unsigned int)*fillevents[i].length, 1, fpw);
evts_tot_write += nFill; evts_tot_write += nFill;
}else{ }else{
if ( count < debugCount) printf("----------------- dropped %d\n", count); if ( count < debugCount) printf("----------------- dropped \n");
evts_tot_drop += nFill; evts_tot_drop += nFill;
} }

431
armory/pxi-time-order.c Normal file
View File

@ -0,0 +1,431 @@
/*********************************************************/
/* 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;
}

16
armory/to2root.cpp
View File

@ -38,6 +38,7 @@
unsigned long long int dataCount=0; unsigned long long int dataCount=0;
unsigned long long int pileUpCount=0; unsigned long long int pileUpCount=0;
unsigned long long int evtCount=0; unsigned long long int evtCount=0;
int traceLength = 280;
/////////////////////////////////// ///////////////////////////////////
// START OF MAIN FUNCTION // // START OF MAIN FUNCTION //
@ -86,9 +87,11 @@ int main(int argc, char **argv) {
unsigned long long e_t[MAX_ID] = {0}; unsigned long long e_t[MAX_ID] = {0};
bool pileup[MAX_ID] = {0}; bool pileup[MAX_ID] = {0};
int qdc[MAX_ID][8] = {0}; int qdc[MAX_ID][8] = {0};
int cfd[MAX_ID] = {0}; //Added by jake to pull out the cfd
int multiCry = 0 ; /// this is total multiplicity for all crystal int multiCry = 0 ; /// this is total multiplicity for all crystal
int runID = 0; // date-run-fileID, Dec15-02-001 = 1502001 int runID = 0; // date-run-fileID, Dec15-02-001 = 1502001
int multiGagg = 0; int multiGagg = 0;
int trace[MAX_ID][traceLength] = {0}; // added by Jake to pull out traces?
//unsigned short pileup[MAXMULTI]; //unsigned short pileup[MAXMULTI];
tree->Branch("evID", &evID, "event_ID/l"); tree->Branch("evID", &evID, "event_ID/l");
@ -98,9 +101,11 @@ int main(int argc, char **argv) {
tree->Branch("e_t", e_t, "e_timestamp[multi]/l"); tree->Branch("e_t", e_t, "e_timestamp[multi]/l");
tree->Branch("pileup", pileup, "pileup[multi]/O"); tree->Branch("pileup", pileup, "pileup[multi]/O");
tree->Branch("qdc", qdc, "qdc[multi][8]/I"); tree->Branch("qdc", qdc, "qdc[multi][8]/I");
tree->Branch("cfd", cfd, "cfd[multi]/I");//added by jake to hold the cfd info
tree->Branch("multiCry", &multiCry, "multiplicity_crystal/I"); tree->Branch("multiCry", &multiCry, "multiplicity_crystal/I");
tree->Branch("multiGagg", &multiGagg, "multiplicity_GAGG/I"); tree->Branch("multiGagg", &multiGagg, "multiplicity_GAGG/I");
tree->Branch("runID", &runID, "runID/I"); tree->Branch("runID", &runID, "runID/I");
tree->Branch("trace", trace, Form("trace[multi][%d]/I",traceLength));
int countGP = 0; //gamma-particle coincident int countGP = 0; //gamma-particle coincident
double totalDataSize = 0; double totalDataSize = 0;
@ -166,8 +171,14 @@ int main(int argc, char **argv) {
id[multi] = mapping[haha]; id[multi] = mapping[haha];
e[multi] = data->energy; e[multi] = data->energy;
e_t[multi] = data->time; e_t[multi] = data->time;
//cfdt[multi] = data->cfd;
pileup[multi] = data->pileup; pileup[multi] = data->pileup;
for( int i = 0; i < 8; i++) qdc[multi][i] = data->QDCsum[i]; for( int i = 0; i < 8; i++) qdc[multi][i] = data->QDCsum[i];
cfd[multi] = data->cfd;
for (int k = 0; k < 280; k++) {
trace[multi][k] = (k < data->trace_length) ? data->trace[k] : 0;
}
multi++ ; multi++ ;
if( id[multi] < 100 ) multiCry ++; if( id[multi] < 100 ) multiCry ++;
if( id[multi] >= 200 ) multiGagg ++; if( id[multi] >= 200 ) multiGagg ++;
@ -180,6 +191,11 @@ int main(int argc, char **argv) {
e_t[multi] = data->time; e_t[multi] = data->time;
pileup[multi] = data->pileup; pileup[multi] = data->pileup;
for( int i = 0; i < 8; i++) qdc[multi][i] = data->QDCsum[i]; for( int i = 0; i < 8; i++) qdc[multi][i] = data->QDCsum[i];
cfd[multi] = data->cfd;
for (int k = 0; k < 280; k++) {
trace[multi][k] = (k < data->trace_length) ? data->trace[k] : 0;
}
multi++ ; multi++ ;
if( id[multi] < 100 ) multiCry ++; if( id[multi] < 100 ) multiCry ++;
if( id[multi] >= 200 ) multiGagg ++; if( id[multi] >= 200 ) multiGagg ++;

375
armory/xia2root.cpp Normal file
View File

@ -0,0 +1,375 @@
#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;
}

38
mapping.h
View File

@ -1,8 +1,11 @@
/************************************ /************************************
Clover : 0 - 99 Clover : 0 - 99
BGO : 100 - 199 BGO : 100 - 199
GAGG : 200 - 299 GAGG A : 200 - 299
ZERO DEGREE : 300 - 399 GAGG B : 300 - 399
ZERO DEGREE : 400 - 499
BEGE : 44 - 45
LaBr3 : 500 - 599
* *********************************/ * *********************************/
#ifndef MAPPING #ifndef MAPPING
@ -16,20 +19,27 @@ ZERO DEGREE : 300 - 399
#define NGAGG 26 #define NGAGG 26
#define NZEROGAGG 2 ///NZERO is used #define NZEROGAGG 2 ///NZERO is used
int mapping[176] ={ //Help for people looking at this mapping:
//rows are pixie slots, colomns correspond to channels. Key for the first row:
//clover1_blue, clover1_black, clover1_green, clover1_red BGO1, clover2_blue, clover2_black, clover2_green, clover2_red, BGO2, clover3_blue, clover3_black, clover3_green, clover3_red, BGO3, EMPTY_CHANNEL
//This also does not take the missing detectors into account, it assumes a detector in each channel.
int mapping[208] ={
//***************** <-- load indicator for EventBuidler //***************** <-- load indicator for EventBuidler
//-0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 //-0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
0, 1, 2, 3, 100, 4, 5, 6, 7, 101, 8, 9, 10, 11, 102, -1, //mod-0 0, 1, 2, 3, 100, 4, 5, 6, 7, 101, 8, 9, 10, 11, 102, 44, //slot 2
12, 13, 14, 15, 103, 16, 17, 18, 19, 104, 20, 21, 22, 23, 105, -1, //mod-1 12, 13, 14, 15, 103, 16, 17, 18, 19, 104, 20, 21, 22, 23, 105, 45, //slot 3
24, 25, 26, 27, 106, 28, 29, 30, 31, 107, 32, 33, 34, 35, 108, -1, //mod-2 24, 25, 26, 27, 106, 28, 29, 30, 31, 107, 32, 33, 34, 35, 108, -1, //slot 4
36, 37, 38, 39, 109, 40, 41, 42, 43, 110, 300, 301, 200, 201, 202, 203, //mod-3 36, 37, 38, 39, 109, 40, 41, 42, 43, 110, 400, 401, -1, -1, -1, -1, //slot 5
204, 205, 206, 207, 208, 209, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, //mod-4 200, 300, 201, 301, 202, 302, 203, 303, 204, 304, 205, 305, 206, 306, 207, 307, //slot 6 Trinity ring 1 (A and B channels alternating) crystal 1 - 8
210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, //mod-5, Ring 4A 208, 308, 209, 309, 210, 310, 211, 311, 212, 312, 213, 313, 214, 314, 215, 315, //slot 7 Ring 2, A and B alternating crystal 1 - 8
260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, //mod-6, Ring 4B 216, 316, 217, 317, 218, 318, 219, 319, 220, 320, 221, 321, 222, 322, 223, 323, //slot 8 Ring 2 first 4 channels ( crystals 9 - 10 A and B ) the rest is Ring 3 crystals 1 - 6 alternating A and B
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //mod-7 224, 324, 225, 325, 226, 326, 227, 327, 228, 328, 229, 329, 230, 330, 231, 331, //slot 9 Ring 3 crystal 7 - 14 A and B alternating
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //mod-8 232, 332, 233, 333, 234, 334, 235, 335, 236, 336, 237, 337, 238, 338, 239, 339, //slot 10 Ring 4 crystals 1 - 8
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //mod-9 240, 340, 241, 341, 242, 342, 243, 343, 244, 344, 245, 345, 246, 346, 247, 347, //slot 11 Ring 4 crystals 9 - 16
-1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, //mod-10 248, 348, 249, 349, 250, 350, 251, 351, 252, 352, 253, 353, 254, 354, 255, 355, //slot 12 Ring 5 crystals 1 - 8
256, 356, 257, 357, 258, 358, 259, 359, 260, 360, 261, 361, 262, 362, 263, 363, //slot 13 Ring 5 crystals 9 -16
500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, //slot 14 LaBr3
//**************** <-- end of mapping indicator EventBuidler //**************** <-- end of mapping indicator EventBuidler
}; };