basic of the pixie2root.cpp
This commit is contained in:
parent
66958a6d2d
commit
b922b4474a
8
.gitignore
vendored
8
.gitignore
vendored
|
|
@ -1,6 +1,14 @@
|
|||
*.root
|
||||
*.evt
|
||||
*.evt.to
|
||||
*.ev2
|
||||
*.log
|
||||
|
||||
xia2root
|
||||
xia2ev2*
|
||||
pixie2root
|
||||
scan
|
||||
|
||||
*.so
|
||||
*.d
|
||||
*.pcm
|
||||
|
|
|
|||
189
Analyzer.C
Normal file
189
Analyzer.C
Normal file
|
|
@ -0,0 +1,189 @@
|
|||
#define Analyzer_cxx
|
||||
|
||||
#include "Analyzer.h"
|
||||
#include <TH2.h>
|
||||
#include <TStyle.h>
|
||||
#include <TH1.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TMath.h>
|
||||
#include <vector>
|
||||
#include <TStopwatch.h>
|
||||
|
||||
//############################################ User setting
|
||||
|
||||
int rawEnergyRange[2] = {100, 2000};
|
||||
|
||||
double BGO_threshold = 100;
|
||||
|
||||
//############################################ end of user setting
|
||||
|
||||
ULong64_t NumEntries = 0;
|
||||
ULong64_t ProcessedEntries = 0;
|
||||
Float_t Frac = 0.1; ///Progress bar
|
||||
TStopwatch StpWatch;
|
||||
|
||||
//############################################ histogram declaration
|
||||
|
||||
TH2F * heVID;
|
||||
TH1F * he[NCLOVER];
|
||||
|
||||
TH2F * hgg[NCLOVER][NCLOVER];
|
||||
|
||||
TH2F * hcoin;
|
||||
|
||||
///----- after calibration
|
||||
TH2F * heCalVID;
|
||||
TH1F * heCal[NCLOVER]; // BGO veto
|
||||
|
||||
|
||||
void Analyzer::Begin(TTree * tree){
|
||||
|
||||
TString option = GetOption();
|
||||
|
||||
NumEntries = tree->GetEntries();
|
||||
|
||||
printf("======================== histogram declaration\n");
|
||||
heVID = new TH2F("heVID", "e vs ID; det ID; e [ch]", NCLOVER, 0, NCLOVER, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
|
||||
heCalVID = new TH2F("heCalVID", Form("eCal vs ID (BGO veto > %.1f); det ID; e [ch]", BGO_threshold), NCLOVER, 0, NCLOVER, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
|
||||
for( int i = 0; i < NCLOVER; i ++){
|
||||
he[i] = new TH1F( Form("he%02d", i), Form("e -%02d", i), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
|
||||
heCal[i] = new TH1F(Form("heCal%02d", i), Form("e -%02d (BGO veto > %.1f)", i, BGO_threshold), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
|
||||
}
|
||||
|
||||
for( int i = 0; i < NCLOVER; i++){
|
||||
for( int j = i; j < NCLOVER; j++){
|
||||
hgg[i][j] = new TH2F(Form("hg%02dg%02d", i, j), Form(" e%02d vs e%02d; e%02d; e%02d", i, j, i, j),
|
||||
rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1],
|
||||
rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
|
||||
}
|
||||
}
|
||||
|
||||
hcoin = new TH2F("hcoin", "detector coin.; det ID; det ID", NCLOVER, 0, NCLOVER, NCLOVER, 0 , NCLOVER);
|
||||
|
||||
printf("======================== End of histograms Declaration\n");
|
||||
StpWatch.Start();
|
||||
|
||||
}
|
||||
|
||||
Bool_t Analyzer::Process(Long64_t entry){
|
||||
|
||||
ProcessedEntries++;
|
||||
|
||||
/*********** Progress Bar ******************************************/
|
||||
if (ProcessedEntries>NumEntries*Frac-1) {
|
||||
TString msg; msg.Form("%llu", NumEntries/1000);
|
||||
int len = msg.Sizeof();
|
||||
printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n",
|
||||
Frac*100, len, ProcessedEntries/1000,NumEntries/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac);
|
||||
StpWatch.Start(kFALSE);
|
||||
Frac+=0.1;
|
||||
}
|
||||
|
||||
b_energy->GetEntry(entry);
|
||||
b_time->GetEntry(entry);
|
||||
b_pileup->GetEntry(entry);
|
||||
b_bgo->GetEntry(entry);
|
||||
b_other->GetEntry(entry);
|
||||
b_multiplicity->GetEntry(entry);
|
||||
|
||||
if( multi == 0 ) return kTRUE;
|
||||
|
||||
for( int detID = 0; detID < NCLOVER ; detID ++){
|
||||
|
||||
//======== baics gate when no energy or pileup
|
||||
if( TMath::IsNaN(e[detID])) continue;
|
||||
//if( pileup[detID] == 1 ) continue;
|
||||
|
||||
//======== Fill raw data
|
||||
heVID->Fill( detID, e[detID]);
|
||||
he[detID]->Fill(e[detID]);
|
||||
|
||||
for( int detJ = detID; detJ < NCLOVER; detJ++) {
|
||||
if( TMath::IsNaN(e[detJ])) continue;
|
||||
hgg[detID][detJ]->Fill(e[detID], e[detJ]); // x then y
|
||||
hcoin->Fill(detID, detJ);
|
||||
}
|
||||
|
||||
//======== BGO veto
|
||||
for( int kk = 0; kk < NBGO ; kk++){
|
||||
if( TMath::IsNaN(bgo[kk]) ) continue;
|
||||
if( bgo[kk] > BGO_threshold ) {
|
||||
return kTRUE;
|
||||
}
|
||||
}
|
||||
|
||||
//========= apply correction
|
||||
double eCal = e[detID];
|
||||
|
||||
heCalVID->Fill( detID, eCal);
|
||||
heCal[detID]->Fill(eCal);
|
||||
|
||||
}
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
|
||||
void listDraws(void) {
|
||||
printf("------------------- List of Plots -------------------\n");
|
||||
printf(" rawID() - Raw e vs ID\n");
|
||||
printf(" drawE() - Raw e for all %d detectors\n", NCLOVER);
|
||||
printf("-----------------------------------------------------\n");
|
||||
}
|
||||
|
||||
void rawID(){
|
||||
TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID");
|
||||
if( cRawID == NULL ) cRawID = new TCanvas("cRawID", "raw ID", 1000, 800);
|
||||
cRawID->cd(1)->SetGrid();
|
||||
heVID->Draw("colz");
|
||||
}
|
||||
|
||||
void drawE(bool isLogy = false, bool cali = false){
|
||||
|
||||
TCanvas *cRawE = (TCanvas *) gROOT->FindObjectAny("cRawE");
|
||||
if( cRawE == NULL ) cRawE = new TCanvas("cRawE", "raw e", 800, 1500);
|
||||
cRawE->Clear();cRawE->Divide(4,9);
|
||||
for (Int_t i = 0; i < 36; i++) {
|
||||
cRawE->cd(i+1);
|
||||
cRawE->cd(i+1)->SetGrid();
|
||||
if( isLogy ) cRawE->cd(i+1)->SetLogy();
|
||||
if( cali ) {
|
||||
heCal[i]->Draw("");
|
||||
}else{
|
||||
he[i]->Draw("");
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
void Analyzer::Terminate(){
|
||||
|
||||
printf("============================== finishing.\n");
|
||||
gROOT->cd();
|
||||
|
||||
int canvasXY[2] = {1200 , 1600} ;// x, y
|
||||
int canvasDiv[2] = {1,2};
|
||||
TCanvas *cCanvas = new TCanvas("cCanvas", "" ,canvasXY[0],canvasXY[1]);
|
||||
cCanvas->Modified(); cCanvas->Update();
|
||||
cCanvas->cd(); cCanvas->Divide(canvasDiv[0],canvasDiv[1]);
|
||||
|
||||
gStyle->SetOptStat("neiou");
|
||||
|
||||
cCanvas->cd(1);
|
||||
cCanvas->cd(1)->SetLogz(1);
|
||||
heVID->Draw("colz");
|
||||
|
||||
cCanvas->cd(2);
|
||||
cCanvas->cd(2)->SetLogz(1);
|
||||
heCalVID->Draw("colz");
|
||||
|
||||
listDraws();
|
||||
|
||||
gROOT->ProcessLine(".L AutoFit.C");
|
||||
printf("=============== loaded AutoFit.C\n");
|
||||
|
||||
|
||||
}
|
||||
115
Analyzer.h
Normal file
115
Analyzer.h
Normal file
|
|
@ -0,0 +1,115 @@
|
|||
//////////////////////////////////////////////////////////
|
||||
// This class has been automatically generated on
|
||||
// Mon Dec 13 18:37:46 2021 by ROOT version 6.24/06
|
||||
// from TTree tree/tree
|
||||
// found on file: efEu152b-000.root
|
||||
//////////////////////////////////////////////////////////
|
||||
|
||||
#ifndef Analyzer_h
|
||||
#define Analyzer_h
|
||||
|
||||
#include <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
|
||||
#include "mapping.h"
|
||||
|
||||
// Header file for the classes stored in the TTree if any.
|
||||
|
||||
class Analyzer : public TSelector {
|
||||
public :
|
||||
TTree *fChain; //!pointer to the analyzed TTree or TChain
|
||||
|
||||
// Fixed size dimensions of array or collections stored in the TTree if any.
|
||||
|
||||
// Declaration of leaf types
|
||||
ULong64_t evID;
|
||||
Double_t e[NCLOVER];
|
||||
ULong64_t t[NCLOVER];
|
||||
UShort_t p[NCLOVER];
|
||||
Double_t bgo[NBGO];
|
||||
Double_t other[NOTHER];
|
||||
Int_t multi;
|
||||
|
||||
// List of branches
|
||||
TBranch *b_event_ID; //!
|
||||
TBranch *b_energy; //!
|
||||
TBranch *b_time; //!
|
||||
TBranch *b_pileup; //!
|
||||
TBranch *b_bgo; //!
|
||||
TBranch *b_other; //!
|
||||
TBranch *b_multiplicity; //!
|
||||
|
||||
Analyzer(TTree * /*tree*/ =0) : fChain(0) { }
|
||||
virtual ~Analyzer() { }
|
||||
virtual Int_t Version() const { return 2; }
|
||||
virtual void Begin(TTree *tree);
|
||||
virtual void SlaveBegin(TTree *tree);
|
||||
virtual void Init(TTree *tree);
|
||||
virtual Bool_t Notify();
|
||||
virtual Bool_t Process(Long64_t entry);
|
||||
virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; }
|
||||
virtual void SetOption(const char *option) { fOption = option; }
|
||||
virtual void SetObject(TObject *obj) { fObject = obj; }
|
||||
virtual void SetInputList(TList *input) { fInput = input; }
|
||||
virtual TList *GetOutputList() const { return fOutput; }
|
||||
virtual void SlaveTerminate();
|
||||
virtual void Terminate();
|
||||
|
||||
ClassDef(Analyzer,0);
|
||||
};
|
||||
|
||||
#endif
|
||||
|
||||
#ifdef Analyzer_cxx
|
||||
void Analyzer::Init(TTree *tree)
|
||||
{
|
||||
// The Init() function is called when the selector needs to initialize
|
||||
// a new tree or chain. Typically here the branch addresses and branch
|
||||
// pointers of the tree will be set.
|
||||
// It is normally not necessary to make changes to the generated
|
||||
// code, but the routine can be extended by the user if needed.
|
||||
// Init() will be called many times when running on PROOF
|
||||
// (once per file to be processed).
|
||||
|
||||
// Set branch addresses and branch pointers
|
||||
if (!tree) return;
|
||||
fChain = tree;
|
||||
fChain->SetMakeClass(1);
|
||||
|
||||
fChain->SetBranchAddress("evID", &evID, &b_event_ID);
|
||||
fChain->SetBranchAddress("e", e, &b_energy);
|
||||
fChain->SetBranchAddress("t", t, &b_time);
|
||||
fChain->SetBranchAddress("p", p, &b_pileup);
|
||||
fChain->SetBranchAddress("bgo", bgo, &b_bgo);
|
||||
fChain->SetBranchAddress("other", other, &b_other);
|
||||
fChain->SetBranchAddress("multi", &multi, &b_multiplicity);
|
||||
}
|
||||
|
||||
Bool_t Analyzer::Notify()
|
||||
{
|
||||
// The Notify() function is called when a new file is opened. This
|
||||
// can be either for a new TTree in a TChain or when when a new TTree
|
||||
// is started when using PROOF. It is normally not necessary to make changes
|
||||
// to the generated code, but the routine can be extended by the
|
||||
// user if needed. The return value is currently not used.
|
||||
|
||||
return kTRUE;
|
||||
}
|
||||
|
||||
|
||||
|
||||
void Analyzer::SlaveBegin(TTree * /*tree*/){
|
||||
|
||||
TString option = GetOption();
|
||||
|
||||
}
|
||||
|
||||
|
||||
void Analyzer::SlaveTerminate(){
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif // #ifdef Analyzer_cxx
|
||||
11
makefile
11
makefile
|
|
@ -1,10 +1,15 @@
|
|||
CC=g++
|
||||
|
||||
all: xia2root
|
||||
#all: xia2root xia2ev2
|
||||
all: xia2root xia2ev2_nopart pixie2root scan
|
||||
|
||||
xia2root: xia2root.cpp
|
||||
$(CC) xia2root.cpp -o xia2root `root-config --cflags --glibs`
|
||||
|
||||
xia2ev2: xia2ev2_nopart.cpp
|
||||
xia2ev2_nopart: xia2ev2_nopart.cpp
|
||||
$(CC) xia2ev2_nopart.cpp -o xia2ev2_nopart
|
||||
|
||||
pixie2root: pixie2root.cpp
|
||||
$(CC) pixie2root.cpp -o pixie2root `root-config --cflags --glibs`
|
||||
|
||||
scan: scan.c
|
||||
$(CC) scan.c -o scan
|
||||
|
|
|
|||
28
mapping.h
Normal file
28
mapping.h
Normal file
|
|
@ -0,0 +1,28 @@
|
|||
/************************************
|
||||
Clover : 0 - 99
|
||||
BGO : 100 - 199
|
||||
Other : 200 - 299
|
||||
|
||||
* *********************************/
|
||||
|
||||
//==================== mapping
|
||||
|
||||
#define NCLOVER 36
|
||||
#define NBGO 9
|
||||
#define NOTHER 62
|
||||
|
||||
// 0 1 2 3 4 5 6 7 8 9
|
||||
int map[130] = { 0, 1, 2, 3, 100, 4, 5, 6, 7, 101, // 0
|
||||
8, 9, 10, 11, 102, -1, 12, 13, 14, 15, // 10
|
||||
103, 16, 17, 18, 19, 104, -1, -1, -1, -1, // 20
|
||||
-1, -1, 20, 21, 22, 23, 105, 24, 25, 26, // 30
|
||||
27, 106, 28, 29, 30, 31, 107, -1, -1, -1, // 40
|
||||
-1, -1, -1, 32, 33, 34, 35, 108, -1, -1, // 50
|
||||
200, 201, 202, 203, 204, 205, 206, 207, 208, 209, // 60
|
||||
210, 211, 212, 213, 214, 215, 216, 217, 218, 219, // 70
|
||||
220, 221, 222, 223, 224, 225, 226, 227, 228, 229, // 80
|
||||
230, 231, 232, 233, 234, 235, 236, 237, 238, 239, // 90
|
||||
240, 241, 242, 243, 244, 245, 246, 247, 248, 249, // 100
|
||||
250, 251, 252, 253, 254, 255, 256, 257, 258, 259, // 110
|
||||
260, 261, -1, -1, -1, -1, -1, -1, -1, -1}; // 120
|
||||
|
||||
474
pixie2root.cpp
Normal file
474
pixie2root.cpp
Normal file
|
|
@ -0,0 +1,474 @@
|
|||
/**********************************************************/
|
||||
/* PXI SCAN CODE -- J.M. Allmond (ORNL) -- July 2016 */
|
||||
/* */
|
||||
/* !unpak data from Pixie-16 digitizers, event build, */
|
||||
/* !and create detctors and user defined spectra */
|
||||
/* */
|
||||
/* gcc -o pxi-scan pxi-scan.c */
|
||||
/* ./pxi-scan -op datafile calibrationfile mapfile */
|
||||
/* */
|
||||
/* ..... calibration file optional */
|
||||
/* ..... map file optional */
|
||||
/* ..... u for update spectra */
|
||||
/* ..... o for overwrite spectra */
|
||||
/* ..... p for print realtime stats */
|
||||
/**********************************************************/
|
||||
|
||||
#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"
|
||||
|
||||
#define PRINT_CAL 1
|
||||
#define PRINT_MAP 1
|
||||
|
||||
#define RAND ((float) rand() / ((unsigned int) RAND_MAX + 1)) // random number in interval (0,1)
|
||||
#define TRUE 1
|
||||
#define FALSE 0
|
||||
|
||||
#define LINE_LENGTH 120
|
||||
|
||||
#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 ; 2004 --> 40 micro-second trace + 4 word header
|
||||
|
||||
#define EVENT_BUILD_TIME 109 // 100 = 1 micro-second ; should be < L + G ~ 5.04 us (note 0.08 us scale factor in set file)
|
||||
|
||||
#define RAWE_REBIN_FACTOR 2.0 // Rebin 32k pixie16 spectra to something smaller to fit better into 8k.
|
||||
|
||||
#include "mapping.h"
|
||||
|
||||
/////////////////////
|
||||
// RAW EVENT TYPES //
|
||||
/////////////////////
|
||||
struct subevent
|
||||
{
|
||||
int chn;
|
||||
int sln;
|
||||
int crn;
|
||||
int id;
|
||||
int hlen;
|
||||
int elen;
|
||||
int trlen; //number of samples
|
||||
int trwlen; //number of words (two samples per word)
|
||||
int fcode; //pileup flag
|
||||
long long int time;
|
||||
int ctime;
|
||||
int ctimef;
|
||||
int energy;
|
||||
int extra;
|
||||
short tr[4096];
|
||||
int esum[4];
|
||||
int qsum[8];
|
||||
};
|
||||
struct subevent subevt[MAX_ID]={0};
|
||||
int sevtmult=0;
|
||||
unsigned long long int sevtcount=0;
|
||||
unsigned long long int pileupcount=0;
|
||||
unsigned long long int evtcount=0;
|
||||
|
||||
int mult[1][4096]={0};
|
||||
|
||||
int tdifid[MAX_ID][8192]={0};
|
||||
|
||||
///////////////////////
|
||||
// Write 2-byte data //
|
||||
///////////////////////
|
||||
void write_data2(char *filename, short *data, int xdim, int ydim, int overwrite) { //2byte per channel Write / Add to previous
|
||||
|
||||
FILE *FP;
|
||||
int i;
|
||||
short *previous;
|
||||
|
||||
if(!overwrite) {
|
||||
//allocate memory for 1d-array for reading in rows of 2d Radware matrix
|
||||
if ( ( previous = (short *)malloc(xdim * ydim * sizeof(short)) ) == NULL ) {
|
||||
printf("\nError, memory not allocated.\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
//open previous spectra file
|
||||
if( (FP=fopen(filename, "r")) != NULL ){
|
||||
fread(previous, sizeof(short)*xdim*ydim, 1, FP);
|
||||
fclose(FP);
|
||||
//update spectra
|
||||
for (i=0; i<xdim*ydim; i++) {
|
||||
if(previous[i] < (powf(2,sizeof(short)*8.0)-2))
|
||||
data[i] = data[i] + previous[i];
|
||||
}
|
||||
}
|
||||
else{
|
||||
printf("%s did not previously exist, creating ...\n", filename);
|
||||
}
|
||||
|
||||
//Deallocate previous data
|
||||
free(previous);
|
||||
}
|
||||
|
||||
FP=fopen(filename, "w");
|
||||
fwrite(data, sizeof(short)*xdim, ydim, FP);
|
||||
fclose(FP);
|
||||
}
|
||||
|
||||
|
||||
///////////////////////
|
||||
// Write 4-byte data //
|
||||
///////////////////////
|
||||
void write_data4(char *filename, int *data, int xdim, int ydim, int overwrite) { //4byte per channel Write / Add to previous
|
||||
|
||||
FILE *FP;
|
||||
int i;
|
||||
int *previous;
|
||||
|
||||
if(!overwrite) {
|
||||
//allocate memory for 1d-array for reading in rows of 2d Radware matrix
|
||||
if ( ( previous = (int *)malloc(xdim * ydim * sizeof(int)) ) == NULL ) {
|
||||
printf("\nError, memory not allocated.\n");
|
||||
exit(1);
|
||||
}
|
||||
|
||||
//open previous spectra file
|
||||
if( (FP=fopen(filename, "r")) != NULL ){
|
||||
fread(previous, sizeof(int)*xdim*ydim, 1, FP);
|
||||
fclose(FP);
|
||||
//update spectra
|
||||
for (i=0; i<xdim*ydim; i++) {
|
||||
if(previous[i] < (powf(2,sizeof(int)*8.0)-2))
|
||||
data[i] = data[i] + previous[i];
|
||||
}
|
||||
}
|
||||
else{
|
||||
printf("%s did not previously exist, creating ...\n", filename);
|
||||
}
|
||||
|
||||
//Deallocate previous data
|
||||
free(previous);
|
||||
}
|
||||
|
||||
FP=fopen(filename, "w");
|
||||
fwrite(data, sizeof(int)*xdim, ydim, FP);
|
||||
fclose(FP);
|
||||
}
|
||||
|
||||
|
||||
///////////////////////////////////
|
||||
// START OF MAIN FUNCTION //
|
||||
///////////////////////////////////
|
||||
int main(int argc, char **argv) {
|
||||
|
||||
int i=0, j=0, k=0;
|
||||
float tempf=0;
|
||||
int max1=0, min1=0;
|
||||
int max2=0, min2=0;
|
||||
int maxid1=-1, minid1=-1;
|
||||
int maxid2=-1, minid2=-1;
|
||||
div_t e_div;
|
||||
lldiv_t lle_div;
|
||||
|
||||
int overwrite = 1;
|
||||
|
||||
double etrace0,etrace1,btrace0,btrace1;
|
||||
double ptrace0,ptrace1,ttrace0,ttrace1,tautrace0,tautrace1;
|
||||
int dbcount = 0;
|
||||
long long int strace[500];
|
||||
memset(strace, 0, sizeof(strace));
|
||||
|
||||
//temp buffer for each sub event
|
||||
unsigned int sub[MAX_SUB_LENGTH];
|
||||
memset(sub, 0, sizeof(sub));
|
||||
|
||||
//Reference time and difference for event building
|
||||
long long int etime, tdif, idtime[MAX_ID]={0}, temptime;
|
||||
|
||||
// Check that the corrent number of arguments were provided.
|
||||
if (argc != 2 && argc != 3 ) {
|
||||
printf("Incorrect number of arguments:\n");
|
||||
printf("%s datafile <outFile>\n", argv[0]);
|
||||
return 1;
|
||||
}
|
||||
|
||||
printf("=====================================\n");
|
||||
printf("=== evt --> root ===\n");
|
||||
printf("=====================================\n");
|
||||
|
||||
//CERN ROOT things
|
||||
TString inFileName = argv[1];
|
||||
TString outFileName = inFileName;
|
||||
|
||||
if( argc >= 3 ){
|
||||
outFileName = argv[2];
|
||||
}else{
|
||||
outFileName.Remove(inFileName.First('.'));
|
||||
outFileName.Append(".root");
|
||||
}
|
||||
printf(" in file : %s \n", inFileName.Data());
|
||||
printf(" our file : %s \n", outFileName.Data());
|
||||
|
||||
printf(" number of detector channal: %d \n", MAX_ID);
|
||||
|
||||
TFile * outRootFile = new TFile(outFileName, "recreate");
|
||||
outRootFile->cd();
|
||||
TTree * tree = new TTree("tree", "tree");
|
||||
|
||||
unsigned long long evID = -1;
|
||||
double energy[NCLOVER];
|
||||
unsigned long long etimestamp[NCLOVER];
|
||||
double bgo[NBGO];
|
||||
double other[NOTHER];
|
||||
unsigned short pileup[NCLOVER];
|
||||
|
||||
//const int maxMulti = 40;
|
||||
//double energy[maxMulti];
|
||||
//unsigned timestamp[maxMulti];
|
||||
//short detID[maxMulti];
|
||||
int multi;
|
||||
|
||||
tree->Branch("evID", &evID, "event_ID/l");
|
||||
///tree->Branch("detID", detID, Form("det ID[%d]/B", NCLOVER));
|
||||
tree->Branch("e", energy, Form("energy[%d]/D", NCLOVER));
|
||||
tree->Branch("t", etimestamp, Form("energy_time_stamp[%d]/l", NCLOVER));
|
||||
tree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCLOVER));
|
||||
|
||||
tree->Branch("bgo", bgo, Form("BGO_energy[%d]/D", NBGO));
|
||||
tree->Branch("other", other, Form("other_energy[%d]/D", NOTHER));
|
||||
|
||||
tree->Branch("multi", &multi, "multiplicity/I");
|
||||
|
||||
//open list-mode data file from PXI digitizer
|
||||
FILE *fpr;
|
||||
long int fprsize,fprpos;
|
||||
if ((fpr = fopen(argv[1], "r")) == NULL) {
|
||||
fprintf(stderr, "Error, cannot open input file %s\n", argv[2]);
|
||||
return 1;
|
||||
}
|
||||
|
||||
//get file size
|
||||
fseek(fpr, 0L, SEEK_END);
|
||||
fprsize = ftell(fpr);
|
||||
rewind(fpr);
|
||||
|
||||
TBenchmark gClock;
|
||||
gClock.Reset();
|
||||
gClock.Start("timer");
|
||||
|
||||
/////////////////////
|
||||
// MAIN WHILE LOOP //
|
||||
/////////////////////
|
||||
while (1) { //main while loop
|
||||
|
||||
/////////////////////////////////
|
||||
// UNPACK DATA AND EVENT BUILD //
|
||||
/////////////////////////////////
|
||||
|
||||
//CERN data clear
|
||||
for( int haha = 0; haha < NCLOVER; haha++){
|
||||
energy[haha] = TMath::QuietNaN();
|
||||
etimestamp[haha] = 0;
|
||||
pileup[haha] = 0;
|
||||
}
|
||||
for( int haha = 0; haha < NBGO; haha++) bgo[haha] = TMath::QuietNaN();
|
||||
for( int haha = 0; haha < NOTHER; haha++) other[haha] = TMath::QuietNaN();
|
||||
multi = 0;
|
||||
evID++;
|
||||
|
||||
etime=-1; tdif=-1; sevtmult=0;
|
||||
//memset(&subevt, 0, sizeof(subevt)); //not needed since everything is redefined (except maybe trace on pileup evts)
|
||||
while (1) { //get subevents and event build for one "event"
|
||||
|
||||
// memset(&subevt[sevtmult], 0, sizeof(subevt[sevtmult])); //not needed since everything is redefined (except maybe trace on pileup evts)
|
||||
|
||||
//read 4-byte header
|
||||
if (fread(sub, sizeof(int)*HEADER_LENGTH, 1, fpr) != 1) break;
|
||||
subevt[sevtmult].chn = sub[0] & 0xF; /// channel in digitizer
|
||||
subevt[sevtmult].sln = (sub[0] & 0xF0) >> 4; /// digitizer ID
|
||||
subevt[sevtmult].crn = (sub[0] & 0xF00) >> 8; /// crate
|
||||
subevt[sevtmult].id = subevt[sevtmult].crn*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (subevt[sevtmult].sln-BOARD_START)*MAX_CHANNELS_PER_BOARD + subevt[sevtmult].chn;
|
||||
subevt[sevtmult].hlen = (sub[0] & 0x1F000) >> 12;
|
||||
subevt[sevtmult].elen = (sub[0] & 0x7FFE0000) >> 17;
|
||||
subevt[sevtmult].fcode = (sub[0] & 0x80000000) >> 31;
|
||||
subevt[sevtmult].time = ( (long long int)(sub[2] & 0xFFFF) << 32) + sub[1];
|
||||
subevt[sevtmult].ctime = (sub[2] & 0x7FFF0000) >> 16;
|
||||
subevt[sevtmult].ctimef = (sub[2] & 0x80000000) >> 31;
|
||||
subevt[sevtmult].energy = (sub[3] & 0xFFFF);
|
||||
subevt[sevtmult].trlen = (sub[3] & 0x7FFF0000) >> 16;
|
||||
subevt[sevtmult].trwlen = subevt[sevtmult].trlen / 2;
|
||||
subevt[sevtmult].extra = (sub[3] & 0x80000000) >> 31;
|
||||
|
||||
//rebin raw trap energy from 32k to ....
|
||||
tempf = (float)subevt[sevtmult].energy/RAWE_REBIN_FACTOR;// + RAND;
|
||||
subevt[sevtmult].energy = (int)tempf;
|
||||
|
||||
//check lengths (sometimes all of the bits for trace length are turned on ...)
|
||||
/* if (subevt[sevtmult].elen - subevt[sevtmult].hlen != subevt[sevtmult].trwlen) {
|
||||
printf("SEVERE ERROR: event, header, and trace length inconsistencies found\n");
|
||||
printf("event length = %d\n", subevt[sevtmult].elen);
|
||||
printf("header length = %d\n", subevt[sevtmult].hlen);
|
||||
printf("trace length = %d\n", subevt[sevtmult].trwlen);
|
||||
printf("Extra = %d\n", subevt[sevtmult].extra);
|
||||
printf("fcode = %d\n", subevt[sevtmult].fcode);
|
||||
//sleep(1);
|
||||
//return 0;
|
||||
} */
|
||||
|
||||
|
||||
//CERN fill tree
|
||||
|
||||
///========== need a mapping, can reduce the array size, speed up.
|
||||
|
||||
int ch = map[subevt[sevtmult].id];
|
||||
if ( 0 <= ch && ch < NCLOVER ){
|
||||
energy[ch] = subevt[sevtmult].energy;
|
||||
etimestamp[ch] = subevt[sevtmult].time;
|
||||
pileup[ch] = subevt[sevtmult].fcode;
|
||||
multi++;
|
||||
}
|
||||
if ( 100 <= ch && ch < 100 + NBGO ){
|
||||
bgo[ch-100] = subevt[sevtmult].energy;
|
||||
}
|
||||
if ( 200 <= ch && ch < 200 + NOTHER){
|
||||
other[ch-200] = subevt[sevtmult].energy;
|
||||
}
|
||||
|
||||
//Set reference time for event building
|
||||
if (etime == -1) {
|
||||
etime = subevt[sevtmult].time;
|
||||
tdif = 0;
|
||||
}
|
||||
else {
|
||||
tdif = subevt[sevtmult].time - etime;
|
||||
if (tdif < 0) {
|
||||
printf("SEVERE ERROR: tdiff < 0, file must be time sorted\n");
|
||||
printf("etime = %lld, time = %lld, and tdif = %lld\n", etime, subevt[sevtmult].time, tdif);
|
||||
return 0;
|
||||
}
|
||||
}
|
||||
|
||||
//Check for end of event, rewind, and break out of while loop
|
||||
if (tdif > EVENT_BUILD_TIME) {
|
||||
fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR); //fwrite/fread is buffered by system ; storing this in local buffer is no faster!
|
||||
break;
|
||||
}
|
||||
|
||||
|
||||
//time between sequential events for a single channel ; useful for determining optimal event building time
|
||||
temptime = (subevt[sevtmult].time - idtime[subevt[sevtmult].id])/100; //rebin to 1 micro-second
|
||||
if ( temptime >= 0 && temptime < 8192) {
|
||||
tdifid[subevt[sevtmult].id][temptime]++;
|
||||
}
|
||||
idtime[subevt[sevtmult].id]=subevt[sevtmult].time; //store time for next subevent of channel
|
||||
|
||||
// total pileups
|
||||
if (subevt[sevtmult].fcode==1) {
|
||||
pileupcount++;
|
||||
}
|
||||
|
||||
|
||||
//more data than just the header; read entire sub event
|
||||
fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR);
|
||||
if (fread(sub, sizeof(int)*subevt[sevtmult].elen, 1, fpr) != 1) break;
|
||||
|
||||
/*
|
||||
//trace
|
||||
k=0;
|
||||
for (i = subevt[sevtmult].hlen; i < subevt[sevtmult].elen; i++) {
|
||||
subevt[sevtmult].tr[i - subevt[sevtmult].hlen + k] = sub[i] & 0x3FFF;
|
||||
subevt[sevtmult].tr[i - subevt[sevtmult].hlen + k + 1] = (sub[i]>>16) & 0x3FFF;
|
||||
k=k+1;
|
||||
}
|
||||
|
||||
// if (subevt[sevtmult].id == 4 && subevt[sevtmult].fcode == 1) DB(subevt[sevtmult].tr);
|
||||
|
||||
//continue if no esum or qsum
|
||||
if (subevt[sevtmult].hlen==HEADER_LENGTH) {
|
||||
sevtmult++;
|
||||
continue;
|
||||
}
|
||||
|
||||
//esum
|
||||
if (subevt[sevtmult].hlen==8 || subevt[sevtmult].hlen==16) {
|
||||
for (i=4; i < 8; i++) {
|
||||
subevt[sevtmult].esum[i-4] = sub[i];
|
||||
}
|
||||
}
|
||||
|
||||
//qsum
|
||||
if (subevt[sevtmult].hlen==12) {
|
||||
for (i=4; i < 12; i++) {
|
||||
subevt[sevtmult].qsum[i-4] = sub[i];
|
||||
}
|
||||
}
|
||||
|
||||
//qsum
|
||||
if (subevt[sevtmult].hlen==16) {
|
||||
for (i=8; i < 16; i++) {
|
||||
subevt[sevtmult].qsum[i-8] = sub[i];
|
||||
}
|
||||
}
|
||||
*/
|
||||
sevtmult++;
|
||||
|
||||
} //end while loop for unpacking sub events and event building for one "event"
|
||||
if (sevtmult==0) break; //end main WHILE LOOP when out of events
|
||||
mult[0][sevtmult]++; //Histogram raw sub event multiplicity
|
||||
sevtcount += sevtmult;
|
||||
evtcount++; //event-built number
|
||||
/////////////////////////////////////
|
||||
// END UNPACK DATA AND EVENT BUILD //
|
||||
/////////////////////////////////////
|
||||
|
||||
//event stats, print status every 10000 events
|
||||
lle_div=lldiv(evtcount,10000);
|
||||
if ( lle_div.rem == 0 ) {
|
||||
fprpos = ftell(fpr);
|
||||
tempf = (float)fprsize/(1024.*1024.*1024.);
|
||||
gClock.Stop("timer");
|
||||
double time = gClock.GetRealTime("timer");
|
||||
gClock.Start("timer");
|
||||
printf("Total SubEvents: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f <mult>)\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[3A\r",
|
||||
sevtcount, (int)((100*pileupcount)/sevtcount), evtcount, (float)sevtcount/(float)evtcount, (100*fprpos/fprsize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.);
|
||||
}
|
||||
|
||||
|
||||
//cern fill tree
|
||||
outRootFile->cd();
|
||||
tree->Fill();
|
||||
|
||||
|
||||
} // end main while loop
|
||||
/////////////////////////
|
||||
// END MAIN WHILE LOOP //
|
||||
/////////////////////////
|
||||
fprpos = ftell(fpr);
|
||||
tempf = (float)fprsize/(1024.*1024.*1024.);
|
||||
printf("Total SubEvents: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f <mult>)\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\n\033[3A\n",
|
||||
sevtcount, (int)((100*pileupcount)/sevtcount), evtcount, (float)sevtcount/(float)evtcount, (100*fprpos/fprsize), tempf);
|
||||
|
||||
//cern save root
|
||||
outRootFile->cd();
|
||||
tree->Write();
|
||||
outRootFile->Close();
|
||||
|
||||
fclose(fpr);
|
||||
|
||||
gClock.Stop("timer");
|
||||
double 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.);
|
||||
|
||||
|
||||
return 0;
|
||||
}
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user