diff --git a/.gitignore b/.gitignore index e561c56..12e1db9 100644 --- a/.gitignore +++ b/.gitignore @@ -25,3 +25,4 @@ Data obsolete test.cpp +test diff --git a/Analyzer.C b/Analyzer.C index 5fb7522..1a19163 100644 --- a/Analyzer.C +++ b/Analyzer.C @@ -36,6 +36,10 @@ TH2F * hcrystalBGO; TH1F * hTDiff; +TH2F * hPID; +TH2F * hPID209; +TH2F * hPID219; + ///----- after calibration and BGO veto TH2F * heCalVID; TH1F * heCal[NCRYSTAL]; @@ -70,6 +74,11 @@ void Analyzer::Begin(TTree * tree){ heCal[i] = new TH1F(Form("heCal%02d", i), Form("eCal -%02d (BGO veto > %.1f); Energy [keV]; count / %d keV", i, BGO_threshold, energyRange[0]), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); } + hPID = new TH2F("hPID", "PID; tail; peak ", 400, -20, 300, 400, -50, 800); + hPID209 = new TH2F("hPID209", "PID detID = 209; tail; peak ", 400, -20, 300, 400, -50, 800); + hPID219 = new TH2F("hPID219", "PID detID = 219; tail; peak ", 400, -20, 300, 400, -50, 800); + + for( int i = 0; i < NCRYSTAL; i++){ for( int j = i+1; j < NCRYSTAL; j++){ //hgg[i][j] = new TH2F(Form("hgg%02d%02d", i, j), Form("e%02d vs e%02d; e%02d; e%02d", i, j, i, j), @@ -112,13 +121,18 @@ Bool_t Analyzer::Process(Long64_t entry){ b_multi->GetEntry(entry); b_multiCry->GetEntry(entry); b_detID->GetEntry(entry); + b_qdc->GetEntry(entry); if( multi == 0 ) return kTRUE; for( int i = 0; i < NCRYSTAL; i++) eCal[i] = TMath::QuietNaN(); ///printf("---------------------------- %d \n", multi); - ///=========== Looping data for the event + + double bgC[2]={0}, peakC[2]={0}, tailC[2] = {0}; + int count = 0 ; + + ///=========== Looping data for the event for( int i = 0; i < multi ; i ++){ int id = detID[i]; @@ -142,6 +156,30 @@ Bool_t Analyzer::Process(Long64_t entry){ } + if( 200 < id && id < 300){ /// GAGG + + } + + if( id == 209 ) { + double bg = (qdc[i][0] + qdc[i][1])/60.; + double peak = qdc[i][3]/20. - bg; + double tail = qdc[i][5]/55. - bg; + hPID209->Fill( tail, peak); + } + if( id == 219 ) { + double bg = (qdc[i][0] + qdc[i][1])/60.; + double peak = qdc[i][3]/20. - bg; + double tail = qdc[i][5]/55. - bg; + hPID219->Fill( tail, peak); + } + + if( count < 2 && (id == 209 || id == 219 )){ + bgC[count] = (qdc[i][0] + qdc[i][1])/60.; + peakC[count] = qdc[i][3]/20. - bgC[count]; + tailC[count] = qdc[i][5]/55. - bgC[count]; + count++; + } + if ( i > 0 ) hTDiff->Fill( e_t[i] - e_t[0]); @@ -182,6 +220,9 @@ Bool_t Analyzer::Process(Long64_t entry){ } + hPID->Fill( (tailC[0]+tailC[1])/2., ( peakC[0] + peakC[1])/2.); + + if ( saveEV2) Save2ev2(); @@ -205,19 +246,20 @@ void Analyzer::Terminate(){ cCanvas->cd(1); cCanvas->cd(1)->SetLogz(1); - heVID->Draw("colz"); + hPID209->Draw("colz"); cCanvas->cd(2); cCanvas->cd(2)->SetLogz(1); - heCalVID->Draw("colz"); + hPID219->Draw("colz"); cCanvas->cd(3); cCanvas->cd(3)->SetLogz(1); - hcrystalBGO->Draw("colz"); + hPID->Draw("colz"); cCanvas->cd(4); cCanvas->cd(4)->SetLogz(1); - hcrystalBGO_G->Draw("colz"); + //gROOT->ProcessLine(".x script.C"); + //hcrystalBGO_G->Draw("colz"); //he[0]->SetLineColor(2); //he[0]->Draw(); //heCal[0]->Draw("same"); diff --git a/Analyzer.h b/Analyzer.h index ea9cf1a..a3a03bd 100644 --- a/Analyzer.h +++ b/Analyzer.h @@ -18,7 +18,9 @@ #include "mapping.h" #include "armory/AnalysisLibrary.h" -// Header file for the classes stored in the TTree if any. +// Header file for the classes stored in the TTree if any.' + +#define MAX_MULTI 200 class Analyzer : public TSelector { public : @@ -28,9 +30,10 @@ public : // Declaration of leaf types ULong64_t evID; - Int_t detID[200]; - Double_t e[200]; - ULong64_t e_t[200]; + Int_t detID[MAX_MULTI]; + Double_t e[MAX_MULTI]; + ULong64_t e_t[MAX_MULTI]; + Int_t qdc[MAX_MULTI][8]; Int_t multi; Int_t multiCry; @@ -39,6 +42,7 @@ public : TBranch *b_detID; //! TBranch *b_energy; //! TBranch *b_time; //! + TBranch *b_qdc; //! TBranch *b_multi; //! TBranch *b_multiCry; //! @@ -99,9 +103,10 @@ void Analyzer::Init(TTree *tree) fChain->SetMakeClass(1); fChain->SetBranchAddress("evID", &evID, &b_event_ID); - fChain->SetBranchAddress("id", detID, &b_detID); + fChain->SetBranchAddress("detID", detID, &b_detID); fChain->SetBranchAddress("e", e, &b_energy); fChain->SetBranchAddress("e_t", e_t, &b_time); + fChain->SetBranchAddress("qdc", qdc, &b_qdc); fChain->SetBranchAddress("multi", &multi, &b_multi); fChain->SetBranchAddress("multiCry", &multiCry, &b_multiCry); diff --git a/armory/DataBlock.h b/armory/DataBlock.h new file mode 100644 index 0000000..59f56ca --- /dev/null +++ b/armory/DataBlock.h @@ -0,0 +1,96 @@ + +#ifndef DATABLOCK_H +#define DATABLOCK_H + +#include +#include +#include +#include +#include + +#include "TSystem.h" +#include "TObject.h" +#include "TFile.h" +#include "TTree.h" +#include "TString.h" + +class DataBlock{ + +public: + UShort_t ch; + UShort_t slot; + UShort_t crate; + UShort_t headerLength; /// headerLength > 4, more data except tarce. + UShort_t eventLength; /// eventLength = headerLength + trace + Bool_t pileup; + ULong64_t time; + UShort_t cfd; + UShort_t energy; + UShort_t trace_length; + Bool_t trace_out_of_range; + + Int_t trailing; + Int_t leading; + Int_t gap; + Int_t baseline; + Int_t QDCsum[8]; + + UShort_t id; + Int_t detID; + ULong64_t eventID; + + UShort_t trace[1024]; + + DataBlock(){ + Clear(); + }; + + ~DataBlock(){}; + + void Clear(){ + ch = 0; + slot = 0; + crate = 0; + headerLength = 0; + eventLength = 0; + pileup = false; + time = 0; + cfd = 0; + energy = 0; + trace_length = 0; + trace_out_of_range = 0; + id = 0; + detID = -1; + eventID = 0; + trailing = 0; + leading = 0; + gap = 0; + baseline = 0; + for( int i = 0; i < 8; i++) QDCsum[i] = -1; + for( int i = 0; i < 1024; i++) trace[i] = 0; + } + + void Print(){ + printf("============== eventID : %llu\n", eventID); + printf("Crate: %d, Slot: %d, Ch: %d | id: %d = detID : %d \n", crate, slot, ch, id, detID); + printf("HeaderLength: %d, Event Length: %d, energy: %d, timeStamp: %llu\n", headerLength, eventLength, energy, time); + printf("trace_length: %d, pile-up:%d\n", trace_length, pileup); + if( headerLength > 4 ){ + if( headerLength > 12 ){ + printf(" trailing : %d\n", trailing); + printf(" leading : %d\n", leading); + printf(" gap : %d\n", gap); + printf(" baseLine : %d\n", baseline); + } + printf(" QDCsum : \n"); + for( int i = 0; i < 8; i++) printf(" %-10d\n", QDCsum[i]); + } + if( eventLength > headerLength ){ + printf(" trace:\n"); + for( int i = 0 ; i < trace_length ; i++)printf("%3d| %-10d\n",i, trace[i]); + } + } + +}; + +#endif diff --git a/armory/MergeEVT.cpp b/armory/MergeEVT.cpp index 42e8677..8ed718e 100644 --- a/armory/MergeEVT.cpp +++ b/armory/MergeEVT.cpp @@ -16,49 +16,8 @@ #define BOARD_START 2 #include "../mapping.h" +#include "../armory/DataBlock.h" -class measurment{ - -public: - UShort_t ch; - UShort_t slot; - UShort_t crate; - UShort_t headerLength; /// headerLength > 4, more data except tarce. - UShort_t eventLength; /// eventLength = headerLength + trace - Bool_t pileup; - ULong64_t time; - UShort_t cfd; - UShort_t energy; - UShort_t trace_length; - Bool_t trace_out_of_range; - - Long64_t timeDiff; - - UShort_t id; - - measurment(){}; - - void Clear(){ - ch = 0; - slot = 0; - crate = 0; - eventLength = 0; - pileup = false; - time = 0; - cfd = 0; - energy = 0; - trace_length = 0; - trace_out_of_range = 0; - timeDiff = 0; - id = 0; - } - - void Print(){ - printf("Crate: %d, Slot: %d, Ch: %d | id: %d \n", crate, slot, ch, id); - printf("HeaderLength: %d, Event Length: %d, energy: %d, timeStamp: %llu\n", headerLength, eventLength, energy, time); - printf("trace_length: %d, pile-up:%d\n", trace_length, pileup); - } -}; //############################################# // main //############################################# @@ -88,7 +47,7 @@ int main(int argn, char **argv) { printf(" out file: %s\n", outFileName.Data()); Long64_t measureID = -1; - measurment data; + DataBlock data; printf("====================================\n"); diff --git a/armory/evt2hist.cpp b/armory/evt2hist.cpp index fe08fa2..e64ae42 100644 --- a/armory/evt2hist.cpp +++ b/armory/evt2hist.cpp @@ -4,6 +4,7 @@ #include #include #include +#include #include "TSystem.h" #include "TObject.h" @@ -22,6 +23,8 @@ #include "../mapping.h" #include "../armory/AnalysisLibrary.h" +#include "../armory/DataBlock.h" + #define MAX_CRATES 2 #define MAX_BOARDS_PER_CRATE 13 #define MAX_CHANNELS_PER_BOARD 16 @@ -32,98 +35,31 @@ // 2) Change to GUI // 4) eventBuilding -class measurment{ +DataBlock dataList[1000]; +int dataCollected = 0; -public: - UShort_t ch; - UShort_t slot; - UShort_t crate; - UShort_t headerLength; /// headerLength > 4, more data except tarce. - UShort_t eventLength; /// eventLength = headerLength + trace - Bool_t pileup; - ULong64_t time; - UShort_t cfd; - UShort_t energy; - UShort_t trace_length; - Bool_t trace_out_of_range; +void BuildEvent(){ - Int_t trailing; - Int_t leading; - Int_t gap; - Int_t baseline; - Int_t QDCsum[8]; + ///==== sort timestamp - UShort_t id; - Int_t detID; - ULong64_t eventID; + ///==== build events - UShort_t trace[1024]; - - measurment(){ - Clear(); - }; - - ~measurment(){}; - - void Clear(){ - ch = 0; - slot = 0; - crate = 0; - headerLength = 0; - eventLength = 0; - pileup = false; - time = 0; - cfd = 0; - energy = 0; - trace_length = 0; - trace_out_of_range = 0; - id = 0; - detID = -1; - eventID = 0; - trailing = 0; - leading = 0; - gap = 0; - baseline = 0; - for( int i = 0; i < 8; i++) QDCsum[i] = -1; - for( int i = 0; i < 1024; i++) trace[i] = 0; - } - - void Print(){ - printf("============== eventID : %llu\n", eventID); - printf("Crate: %d, Slot: %d, Ch: %d | id: %d = detID : %d \n", crate, slot, ch, id, detID); - printf("HeaderLength: %d, Event Length: %d, energy: %d, timeStamp: %llu\n", headerLength, eventLength, energy, time); - printf("trace_length: %d, pile-up:%d\n", trace_length, pileup); - if( headerLength > 4 ){ - if( headerLength > 12 ){ - printf(" trailing : %d\n", trailing); - printf(" leading : %d\n", leading); - printf(" gap : %d\n", gap); - printf(" baseLine : %d\n", baseline); - } - printf(" QDCsum : \n"); - for( int i = 0; i < 8; i++) printf(" %-10d\n", QDCsum[i]); - } - if( eventLength > headerLength ){ - printf(" trace:\n"); - for( int i = 0 ; i < trace_length ; i++)printf("%3d| %-10d\n",i, trace[i]); - } - } - -}; + ///==== output to g-g hist +} //############################################# // main //############################################### int main(int argn, char **argv) { - if (argn < 2 || argn > 6 ) { + if (argn < 2 || argn > 7 ) { printf("Usage :\n"); - printf("%s [evt File] [E corr] [raw E threshold] [Save Hist] [Save Root]\n", argv[0]); + printf("%s [evt File] [E corr] [raw E threshold] [Hist File] [Root File] [display data]\n", argv[0]); printf(" [E corr] : correction file for gamma energy \n"); - printf(" [raw E threshold] : min raw E \n"); - printf(" [save Hist] : 1/0 \n"); - printf(" [save Root] : 1/0 \n"); - + printf(" [raw E threshold] : min raw E, default 10 ch \n"); + printf(" [Hist File] : if provided, saved histograms in root. if value = 1, auto fileName. \n"); + printf(" [Root File] : if provided, saved energy, timestamp, QDC, and trace in to root. if value = 1, auto fileName. \n"); + printf(" [display data] : number of event from the first one to display. default 0. \n"); return 1; } @@ -137,15 +73,31 @@ int main(int argn, char **argv) { eCorr = LoadCorrectionParameters(corrFile); } - int rawEnergyThreshold = 100; + int rawEnergyThreshold = 10; if( argn >= 4 ) rawEnergyThreshold = atoi(argv[3]); - bool isSaveHist = false; ///save gamma hist for calibration - if( argn >= 5 ) isSaveHist = atoi(argv[5]); - - bool isSaveRoot = false; ///save data into root - if( argn >= 6 ) isSaveRoot = atoi(argv[6]); + TString histFileName = ""; ///save gamma hist for calibration + if( argn >= 5 ) histFileName = argv[4]; + if( histFileName == "1" ){ + histFileName = inFileName; + histFileName.Remove(0, inFileName.Last('/')+1); + histFileName.Remove(histFileName.First('.')); + histFileName.Append("_hist.root"); + } + if( histFileName == "0" ) histFileName = ""; + TString rootFileName = ""; ///save data into root + if( argn >= 6 ) rootFileName = argv[5]; + if( rootFileName == "1" ){ + rootFileName = inFileName; + rootFileName.Remove(0, inFileName.Last('/')+1); + rootFileName.Remove(rootFileName.First('.')); + rootFileName.Append("_raw.root"); + } + if( rootFileName == "0" ) rootFileName = ""; + + int maxDataDisplay = 0; + if( argn >= 7 ) maxDataDisplay = atoi(argv[6]); long int inFilePos; TBenchmark gClock; @@ -154,7 +106,7 @@ int main(int argn, char **argv) { ULong64_t measureID = -1; - measurment data; + DataBlock data; printf("====================================\n"); @@ -166,20 +118,27 @@ int main(int argn, char **argv) { printf(" in file: \033[1;31m%s\033[m\n", inFileName.Data()); printf(" Gamma energy correction file : %s\n", corrFile == "" ? "Not provided." : corrFile.Data()); + printf(" raw E threshold : %d ch\n", rawEnergyThreshold); + if( histFileName != "" ) printf(" Save histograms to %s\n", histFileName.Data()); + if( rootFileName != "" ) printf(" Save root to %s\n", rootFileName.Data()); printf("--------------------------------\n"); TFile * fFile = NULL; TTree * tree = NULL; - if( isSaveRoot ){ - fFile = new TFile("temp.root", "RECREATE"); + if( rootFileName != "" ){ + fFile = new TFile(rootFileName, "RECREATE"); tree = new TTree("tree", "tree"); + + tree->Branch("headerLenght", &data.headerLength, "HeaderLength/s"); tree->Branch("detID", &data.detID, "detID/s"); + tree->Branch("id", &data.id, "id/s"); tree->Branch("e", &data.energy, "energy/s"); tree->Branch("e_t", &data.time, "timestamp/l"); + tree->Branch("p", &data.pileup, "pileup/O"); tree->Branch("qdc", data.QDCsum, "QDC_sum[8]/I"); tree->Branch("trace_length", &data.trace_length, "trace_length/s"); - tree->Branch("trace", data.trace, "trace[trace_length]/s"); + tree->Branch("trace", data.trace, "trace[trace_length]/s"); } //================ get file size @@ -227,7 +186,7 @@ int main(int argn, char **argv) { while ( ! feof(inFile) ){ - fread(header, sizeof(header), 1, inFile); + if ( fread(header, sizeof(header), 1, inFile) !=1 ) break; measureID ++; /// see the Pixie-16 user manual, Table4-2 @@ -251,17 +210,18 @@ int main(int argn, char **argv) { if( data.headerLength >= 4 ){ unsigned int extraHeader[data.headerLength-4]; fread(extraHeader, sizeof(extraHeader), 1, inFile); - if( data.headerLength > 12){ + if( data.headerLength == 8 || data.headerLength == 16){ data.trailing = extraHeader[0]; data.leading = extraHeader[1]; data.gap = extraHeader[2]; data.baseline = extraHeader[3]; } - - for( int i = 0; i < 8; i++){ - int startID = 0; - if( data.headerLength > 12) startID = 4; ///the 1st 4 words - data.QDCsum[i] = extraHeader[i+startID]; + if( data.headerLength == 12 || data.headerLength == 16){ + for( int i = 0; i < 8; i++){ + int startID = 0; + if( data.headerLength > 12) startID = 4; ///the 1st 4 words + data.QDCsum[i] = extraHeader[i+startID]; + } } } ///====== read trace @@ -275,11 +235,11 @@ int main(int argn, char **argv) { } } - ///if( measureID < 10 ) { - /// printf("----------------------event Length: %u, fpos: %llu byte (%lld words)\n", data.eventLength, fpos, fpos/4); - /// for(int i = 0; i < 4; i++) printf(" %x\n", header[i]); - /// data.Print(); - ///} + if( measureID < maxDataDisplay ) { + printf("----------------------event Length: %u, fpos: %llu byte (%lld words)\n", data.eventLength, fpos, fpos/4); + for(int i = 0; i < 4; i++) printf(" %x\n", header[i]); + data.Print(); + } //=== jump to next measurement. obsolete, we read the whole block ///if( data.eventLength > 4 ) { @@ -309,12 +269,29 @@ int main(int argn, char **argv) { gTrace->Clear(); gTrace->Set(data.trace_length); gTrace->SetTitle(Form("eventID : %llu, detID: %d", data.eventID, data.detID)); + + if( data.headerLength == 4 ) { + for( int i = 0; i < 8; i++ ) data.QDCsum[i] = 0; + } + for( int i = 0; i < data.trace_length; i++){ gTrace->SetPoint(i, i, data.trace[i]); + + ///if the header don't have ADC, make one + if( data.headerLength < 12 ) { + if( 0 <= i && i < 31 ) data.QDCsum[0] += data.trace[i]; + if( 31 <= i && i < 60 ) data.QDCsum[1] += data.trace[i]; + if( 60 <= i && i < 75 ) data.QDCsum[2] += data.trace[i]; + if( 75 <= i && i < 95 ) data.QDCsum[3] += data.trace[i]; + if( 95 <= i && i < 105 ) data.QDCsum[4] += data.trace[i]; + if( 105 <= i && i < 160 ) data.QDCsum[5] += data.trace[i]; + if( 160 <= i && i < 175 ) data.QDCsum[6] += data.trace[i]; + if( 175 <= i && i < 200 ) data.QDCsum[7] += data.trace[i]; + } } } - if( isSaveRoot ){ + if( rootFileName != "" ){ fFile->cd(); tree->Fill(); } @@ -415,15 +392,15 @@ int main(int argn, char **argv) { printf("\n\n\n============= reached end of file\n"); - if( isSaveHist ) { - printf(" save gamma histograms : \033[1;3mhist.root\033[m\n"); - TFile * fHist = new TFile("hist.root", "RECREATE"); + if( histFileName != "" ) { + printf(" save gamma histograms : \033[1;3m%s\033[m\n", histFileName.Data()); + TFile * fHist = new TFile(histFileName, "RECREATE"); for( int i = 0; i < NCRYSTAL; i++) he[i]->Write("", TObject::kOverwrite); fHist->Close(); } - if( isSaveRoot){ - printf(" save into Root : \033[1;3mtemp.root\033[m\n"); + if( rootFileName != "" ){ + printf(" save into Root : \033[1;3m%s\033[m\n", rootFileName.Data()); fFile->cd(); tree->Write(); fFile->Close(); diff --git a/armory/evtReader.h b/armory/evtReader.h new file mode 100644 index 0000000..f5379ef --- /dev/null +++ b/armory/evtReader.h @@ -0,0 +1,195 @@ +#ifndef EVTREADER_H +#define EVTREADER_H + +#include +#include +#include +#include +#include + +#include "TSystem.h" +#include "TObject.h" +#include "TFile.h" +#include "TTree.h" +#include "TString.h" +#include "TBenchmark.h" + +#include "../mapping.h" +#include "../armory/DataBlock.h" + +#define MAX_CRATES 2 +#define MAX_BOARDS_PER_CRATE 13 +#define MAX_CHANNELS_PER_BOARD 16 +#define BOARD_START 2 + +class evtReader{ + +public: + DataBlock * data; + + Long64_t blockID; + + bool isOpen; + +private: + FILE * inFile; + + long int inFileSize; + long int inFilePos; + bool endOfFile; + + TBenchmark gClock; + +///============================================ Methods +public: + evtReader(TString inFileName){ + + inFile = fopen(inFileName, "r"); + if( inFile == NULL ){ + printf("Cannot read file : %s \n", inFileName.Data()); + + inFileSize = 0; + inFilePos = 0; + + blockID = -1; + endOfFile = false; + isOpen = false; + + }else{ + + fseek(inFile, 0L, SEEK_END); + inFileSize = ftell(inFile); + inFilePos = 0; + rewind(inFile); ///back to the File begining + + data = new DataBlock(); + data->Clear(); + blockID = -1; + + endOfFile = false; + isOpen = true; + + gClock.Reset(); + gClock.Start("timer"); + } + } + + + ~evtReader(){ + fclose(inFile); + delete inFile; + delete data; + } + + void UpdateFileSize(){ + if( inFile == NULL ) return; + fseek(inFile, 0L, SEEK_END); + inFileSize = ftell(inFile); + fseek(inFile, inFilePos, SEEK_SET); + } + + bool isEndOfFile() { + int haha = feof(inFile); + return haha > 0 ? true: false; + } + + long int GetFilePos(){return inFilePos;} + long int GetFileSize(){return inFileSize;} + + int ReadBlock(){ + + if( feof(inFile) ) return -1; + if( endOfFile ) return -1; + + unsigned int header[4]; //read 4 header, unsigned int = 4 byte = 32 bits. + + if ( fread(header, sizeof(header), 1, inFile) != 1 ) { + endOfFile = true; + return -1; + } + blockID ++; + + /// see the Pixie-16 user manual, Table4-2 + data->eventID = blockID; + data->ch = header[0] & 0xF ; + data->slot = (header[0] >> 4) & 0xF; + data->crate = (header[0] >> 8) & 0xF; + data->headerLength = (header[0] >> 12) & 0x1F; + data->eventLength = (header[0] >> 17) & 0x3FFF; + data->pileup = header[0] >> 31 ; + data->time = ((ULong64_t)(header[2] & 0xFFFF) << 32) + header[1]; + data->cfd = header[2] >> 16 ; + data->energy = (header[3] & 0xFFFF ) /2; // I don;t know why it has to "rebin" + data->trace_length = (header[3] >> 16) & 0x7FFF; + data->trace_out_of_range = header[3] >> 31; + + data->id = data->crate*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data->slot-BOARD_START)*MAX_CHANNELS_PER_BOARD + data->ch; + data->detID = mapping[data->id]; + + ///======== read QDCsum + if( data->headerLength >= 4 ){ + unsigned int extraHeader[data->headerLength-4]; + fread(extraHeader, sizeof(extraHeader), 1, inFile); + if( data->headerLength == 8 || data->headerLength == 16){ + data->trailing = extraHeader[0]; + data->leading = extraHeader[1]; + data->gap = extraHeader[2]; + data->baseline = extraHeader[3]; + } + if( data->headerLength == 12 || data->headerLength == 16){ + for( int i = 0; i < 8; i++){ + int startID = 0; + if( data->headerLength > 12) startID = 4; ///the 1st 4 words + data->QDCsum[i] = extraHeader[i+startID]; + } + } + } + ///====== read trace + if( data->eventLength > data->headerLength ){ + unsigned int traceBlock[data->trace_length / 2]; + fread(traceBlock, sizeof(traceBlock), 1, inFile); + + for( int i = 0; i < data->trace_length/2 ; i++){ + data->trace[2*i+0] = traceBlock[i] & 0xFFFF ; + data->trace[2*i+1] = (traceBlock[i] >> 16 ) & 0xFFFF ; + } + + ///make QDC by trace + if( data->headerLength == 8 ) { + for( int i = 0; i < data->trace_length; i++){ + if( 0 <= i && i < 31 ) data->QDCsum[0] += data->trace[i]; + if( 31 <= i && i < 60 ) data->QDCsum[1] += data->trace[i]; + if( 60 <= i && i < 75 ) data->QDCsum[2] += data->trace[i]; + if( 75 <= i && i < 95 ) data->QDCsum[3] += data->trace[i]; + if( 95 <= i && i < 105 ) data->QDCsum[4] += data->trace[i]; + if( 105 <= i && i < 160 ) data->QDCsum[5] += data->trace[i]; + if( 160 <= i && i < 175 ) data->QDCsum[6] += data->trace[i]; + if( 175 <= i && i < 200 ) data->QDCsum[7] += data->trace[i]; + } + } + } + + inFilePos = ftell(inFile); + + return 1; + } + + + void PrintStatus(int id){ + + ///==== event stats, print status every 10000 events + if ( blockID % id == 0 ) { + UpdateFileSize(); + //inFilePos = ftell(inFile); + gClock.Stop("timer"); + double time = gClock.GetRealTime("timer"); + gClock.Start("timer"); + printf("Total measurements: \x1B[32m%llu \x1B[0m\nReading Pos: \x1B[32m %.3f/%.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[A\033[A\r", + blockID, inFilePos/(1024.*1024.*1024.), inFileSize/1024./1024./1024, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + } + + } + +}; + +#endif diff --git a/armory/makefile b/armory/makefile index d3e03af..d24297e 100644 --- a/armory/makefile +++ b/armory/makefile @@ -1,11 +1,6 @@ CC=g++ -#all: xia2root xia2ev2_nopart pixie2root scan pxi-time-order -#all: xia2root xia2ev2_nopart pixie2root scan evt2root evt2hist -#all: xia2root to2root MergeEVT evt2hist pxi-time-order ev22txt EventBuilder -#all: xia2root to2root MergeEVT pxi-time-order ev22txt EventBuilder -#all: to2root MergeEVT ev22txt EventBuilder pxi-time-order -all: to2root evt2hist MergeEVT ev22txt EventBuilder pxi-time-order +all: to2root evt2hist MergeEVT ev22txt EventBuilder pxi-time-order test #this is FSU evt to root xia2root: ../armory/xia2root.cpp @@ -19,11 +14,11 @@ to2root: ../armory/to2root.cpp $(CC) ../armory/to2root.cpp -o to2root `root-config --cflags --glibs` #this is for online root -MergeEVT: ../armory/MergeEVT.cpp +MergeEVT: ../armory/MergeEVT.cpp ../armory/DataBlock.h ../mapping.h $(CC) ../armory/MergeEVT.cpp -o MergeEVT `root-config --cflags --glibs` #this is for online spectrums -evt2hist: ../armory/evt2hist.cpp +evt2hist: ../armory/evt2hist.cpp ../armory/DataBlock.h ../mapping.h $(CC) ../armory/evt2hist.cpp -o evt2hist `root-config --cflags --glibs` pxi-time-order: ../armory/pxi-time-order.c @@ -35,5 +30,8 @@ ev22txt: ../armory/ev22txt.cpp EventBuilder: ../armory/EventBuilder.cpp $(CC) ../armory/EventBuilder.cpp -o EventBuilder `root-config --cflags --glibs` +test: ../armory/test.cpp ../armory/DataBlock.h ../armory/evtReader.h ../mapping.h + $(CC) ../armory/test.cpp -o test `root-config --cflags --glibs` + clean: - -rm xia2root to2root MergeEVT evt2hist pxi-time-order ev22txt EventBuilder + -rm xia2root to2root MergeEVT evt2hist pxi-time-order ev22txt EventBuilder test diff --git a/armory/to2root.cpp b/armory/to2root.cpp index 9e8fa7f..87783ea 100644 --- a/armory/to2root.cpp +++ b/armory/to2root.cpp @@ -121,6 +121,7 @@ int main(int argc, char **argv) { int id[MAX_ID] = {0}; double e[MAX_ID] = {TMath::QuietNaN()}; unsigned long long e_t[MAX_ID] = {0}; + int qdc[MAX_ID][8] = {0}; Int_t multiCry = 0 ; /// this is total multiplicity for all crystal //unsigned short pileup[MAXMULTI]; @@ -128,9 +129,10 @@ int main(int argc, char **argv) { tree->Branch("evID", &evID, "event_ID/l"); tree->Branch("multi", &multi, "multi/I"); - tree->Branch("id", id, "id[multi]/I"); + tree->Branch("detID", id, "detID[multi]/I"); tree->Branch("e", e, "e[multi]/D"); tree->Branch("e_t", e_t, "e_timestamp[multi]/l"); + tree->Branch("qdc", qdc, "qdc[multi][8]/I"); tree->Branch("multiCry", &multiCry, "multiplicity_crystal/I"); //open list-mode data file from PXI digitizer @@ -191,6 +193,49 @@ int main(int argc, char **argv) { tempf = (float)data.e/RAWE_REBIN_FACTOR;// + RAND; data.e = (int)tempf; + if( data.hlen > 4 ){ + unsigned int extraHeader[data.hlen - 4 ]; + fread(extraHeader, sizeof(extraHeader), 1, fpr); + if( data.hlen == 8 || data.hlen == 16){ + for( int i = 0 ; i < 4; i ++) data.esum[i] = extraHeader[i]; + } + + if( data.hlen == 12 || data.hlen == 16){ + for( int i = 0; i < 8; i++){ + int startID = 0; + if( data.hlen > 12) startID = 4; ///the 1st 4 words + data.qsum[i] = extraHeader[i+startID]; + } + } + } + + if( data.hlen < data.elen ){ + unsigned int traceBlock[data.trlen / 2]; + fread(traceBlock, sizeof(traceBlock), 1, fpr); + + for( int i = 0; i < data.trlen/2 ; i++){ + data.tr[2*i+0] = traceBlock[i] & 0xFFFF ; + data.tr[2*i+1] = (traceBlock[i] >> 16 ) & 0xFFFF ; + } + + + if( data.hlen < 12 ){ + for( int i = 0 ; i < 8; i++) data.qsum[i] = 0; + for( int i = 0; i < data.trlen ; i++){ + if( 0 <= i && i < 31 ) data.qsum[0] += data.tr[i]; + if( 31 <= i && i < 60 ) data.qsum[1] += data.tr[i]; + if( 60 <= i && i < 75 ) data.qsum[2] += data.tr[i]; + if( 75 <= i && i < 95 ) data.qsum[3] += data.tr[i]; + if( 95 <= i && i < 105 ) data.qsum[4] += data.tr[i]; + if( 105 <= i && i < 160 ) data.qsum[5] += data.tr[i]; + if( 160 <= i && i < 175 ) data.qsum[6] += data.tr[i]; + if( 175 <= i && i < 200 ) data.qsum[7] += data.tr[i]; + } + } + } + + + //Set reference time for event building if (etime == -1) { etime = data.time; @@ -202,7 +247,7 @@ int main(int argc, char **argv) { //Check for end of event, rewind, and break out of while loop if (tdif > timeWindow) { - fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR); //fwrite/fread is buffered by system ; storing this in local buffer is no faster! + fseek(fpr, -sizeof(int)*data.elen, SEEK_CUR); //fwrite/fread is buffered by system ; storing this in local buffer is no faster! break; }else{ //if within time window, fill array; @@ -210,6 +255,7 @@ int main(int argc, char **argv) { id[multi] = detID; e[multi] = data.e; e_t[multi] = data.time; + for( int i = 0; i < 8; i++) qdc[multi][i] = data.qsum[i]; multi++ ; if( detID < 100 ) multiCry ++; } @@ -220,49 +266,9 @@ int main(int argc, char **argv) { } //more data than just the header; read entire sub event, first rewind, then read data.elen - fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR); + //fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR); //if (fread(sub, sizeof(int)*dataBlock[sevtmult].elen, 1, fpr) != 1) break; - if (fread(sub, sizeof(int)*data.elen, 1, fpr) != 1) break; - - - /** - //trace - k=0; - for (i = dataBlock[sevtmult].hlen; i < dataBlock[sevtmult].elen; i++) { - dataBlock[sevtmult].tr[i - dataBlock[sevtmult].hlen + k] = sub[i] & 0x3FFF; - dataBlock[sevtmult].tr[i - dataBlock[sevtmult].hlen + k + 1] = (sub[i]>>16) & 0x3FFF; - k=k+1; - } - - // if (dataBlock[sevtmult].id == 4 && dataBlock[sevtmult].fcode == 1) DB(dataBlock[sevtmult].tr); - - //continue if no esum or qsum - if (dataBlock[sevtmult].hlen==HEADER_LENGTH) { - sevtmult++; - continue; - } - - //esum - if (dataBlock[sevtmult].hlen==8 || dataBlock[sevtmult].hlen==16) { - for (i=4; i < 8; i++) { - dataBlock[sevtmult].esum[i-4] = sub[i]; - } - } - - //qsum - if (dataBlock[sevtmult].hlen==12) { - for (i=4; i < 12; i++) { - dataBlock[sevtmult].qsum[i-4] = sub[i]; - } - } - - //qsum - if (dataBlock[sevtmult].hlen==16) { - for (i=8; i < 16; i++) { - dataBlock[sevtmult].qsum[i-8] = sub[i]; - } - } - */ + //if (fread(sub, sizeof(int)*data.elen, 1, fpr) != 1) break; } //end while loop for unpacking sub events and event building for one "event" if (multi==0) break; //end main WHILE LOOP when out of events diff --git a/mapping.h b/mapping.h index 11183a8..e9d272c 100644 --- a/mapping.h +++ b/mapping.h @@ -4,6 +4,8 @@ BGO : 100 - 199 Other : 200 - 299 * *********************************/ +#ifndef MAPPING +#define MAPPING //==================== mapping @@ -27,3 +29,4 @@ int mapping[130] ={ 0, 1, 2, 3, 100, 4, 5, 6, 7, 101, // 0 250, 251, -1, -1, -1, -1, -1, -1, -1, -1, // 110 -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // 120 +#endif