From 60af179708d60f3e4f82fdf1caeca15e892da35d Mon Sep 17 00:00:00 2001 From: "Ryan@WorkStation" Date: Thu, 6 Jan 2022 13:18:29 -0500 Subject: [PATCH] complete all basic programs for single evt to _raw.root and root --- .gitignore | 6 +- Analyzer.C | 123 +++++++---- Analyzer.h | 47 ++++- armory/Analyzer_Utili.c | 12 +- armory/EventBuilder.cpp | 203 ++++++++++++++++++ armory/ev22txt.cpp | 81 +++++++ armory/evt2root.cpp | 27 +-- armory/pixie2root.cpp | 454 ---------------------------------------- armory/to2root.cpp | 380 +++++++++++++++++++++++++++++++++ makefile | 12 +- mapping.h | 28 +-- 11 files changed, 836 insertions(+), 537 deletions(-) create mode 100644 armory/EventBuilder.cpp create mode 100644 armory/ev22txt.cpp delete mode 100644 armory/pixie2root.cpp create mode 100644 armory/to2root.cpp diff --git a/.gitignore b/.gitignore index b889b66..9719c91 100644 --- a/.gitignore +++ b/.gitignore @@ -6,11 +6,12 @@ xia2root xia2ev2* -pixie2root +to2root scan pxi-time-order evt2root evt2hist +ev22txt *.so *.d @@ -18,3 +19,6 @@ evt2hist ti74pt7a fsu_run_2021 +data + +test.cpp diff --git a/Analyzer.C b/Analyzer.C index 60df793..322d601 100644 --- a/Analyzer.C +++ b/Analyzer.C @@ -7,25 +7,19 @@ #include #include #include -#include +#include //############################################ User setting -int rawEnergyRange[2] = {500, 6000}; // in ch -int energyRange[3] = {1, 100, 2000}; // keV {resol, min, max} +int rawEnergyRange[2] = {0, 6000}; // in ch +int energyRange[3] = {1, 0, 6000}; // keV {resol, min, max} -double BGO_threshold = 100; // in ch +double BGO_threshold = 0; // in ch -TString e_corr = "correction_e.dat"; +TString e_corr = "";//"correction_e.dat"; //############################################ end of user setting -ULong64_t NumEntries = 0; -ULong64_t ProcessedEntries = 0; -Float_t Frac = 0.1; ///Progress bar -TStopwatch StpWatch; - -vector> eCorr; //############################################ histogram declaration @@ -41,12 +35,19 @@ TH2F * heCalvID; TH1F * heCal[NCRYSTAL]; TH2F * hcoinBGO; +TH2F * hcrystalBGO; + //############################################ BEGIN void Analyzer::Begin(TTree * tree){ TString option = GetOption(); - NumEntries = tree->GetEntries(); + totnumEntry = tree->GetEntries(); + + printf( "=========================================================================== \n"); + printf( "========================== Analysis.C/h ================================ \n"); + printf( "====== total Entry : %lld \n", totnumEntry); + printf( "=========================================================================== \n"); printf("======================== Histograms declaration\n"); @@ -70,16 +71,15 @@ void Analyzer::Begin(TTree * tree){ } hcoin = new TH2F("hcoin", "detector coin.; det ID; det ID", NCRYSTAL, 0, NCRYSTAL, NCRYSTAL, 0 , NCRYSTAL); + hcoinBGO = new TH2F("hcoinBGO", Form("detector coin. (BGO veto > %.1f); det ID; det ID", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, NCRYSTAL, 0 , NCRYSTAL); + hcrystalBGO = new TH2F("hcrystalBGO", Form("crystal vs BGO ; det ID; BGO ID"), NCRYSTAL, 0, NCRYSTAL, NBGO, 0 , NBGO); printf("======================== End of histograms declaration\n"); printf("======================== Load parameters.\n"); eCorr = LoadCorrectionParameters(e_corr); - - StpWatch.Start(); - printf("======================== Start processing....\n"); } //############################################ PROCESS @@ -88,23 +88,30 @@ Bool_t Analyzer::Process(Long64_t entry){ ProcessedEntries++; /*********** Progress Bar ******************************************/ - if (ProcessedEntries>NumEntries*Frac-1) { - TString msg; msg.Form("%llu", NumEntries/1000); + if (ProcessedEntries>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, ProcessedEntries/1000,NumEntries/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); + Frac*100, len, ProcessedEntries/1000,totnumEntry/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_pileup->GetEntry(entry); + b_hit->GetEntry(entry); b_bgo->GetEntry(entry); - //b_other->GetEntry(entry); - //b_multiplicity->GetEntry(entry); + b_bgoTime->GetEntry(entry); + b_other->GetEntry(entry); + b_multi->GetEntry(entry); if( multi == 0 ) return kTRUE; + + int numGatedData=0; + vector gatedID; + gatedID.clear(); + double eCal[NCRYSTAL] ={0.0}; ///=========== Looping Crystals for( int detID = 0; detID < NCRYSTAL ; detID ++){ @@ -117,6 +124,10 @@ Bool_t Analyzer::Process(Long64_t entry){ hevID->Fill( detID, e[detID]); he[detID]->Fill(e[detID]); + for( int kk = 0 ; kk < NBGO; kk++){ + if( bgo[kk] > 0 ) hcrystalBGO->Fill(detID, kk); + } + for( int detJ = detID +1; detJ < NCRYSTAL; detJ++) { if( TMath::IsNaN(e[detJ])) continue; @@ -127,34 +138,73 @@ Bool_t Analyzer::Process(Long64_t entry){ //======== BGO veto for( int kk = 0; kk < NBGO ; kk++){ if( TMath::IsNaN(bgo[kk]) ) continue; - if( bgo[kk] > BGO_threshold ) { + if( bgo[kk] > BGO_threshold && 4*kk <= detID && detID < 4*(kk+1) ) { return kTRUE; } } - - //========= apply correction - double eCal = 0 ; - int order = (int) eCorr[detID].size(); - for( int i = 0; i < order ; i++){ - eCal += eCorr[detID][i] * TMath::Power(e[detID], i); + //----- for ev2 file + if( saveEV2 ){ + numGatedData ++; + gatedID.push_back(detID); } - - heCalvID->Fill( detID, eCal); - heCal[detID]->Fill(eCal); + //========= apply correction + int order = (int) eCorr[detID].size(); + for( int i = 0; i < order ; i++){ + eCal[detID] += eCorr[detID][i] * TMath::Power(e[detID], i); + } + + heCalvID->Fill( detID, eCal[detID]); + heCal[detID]->Fill(eCal[detID]); for( int detJ = detID +1; detJ < NCRYSTAL; detJ++) { if( TMath::IsNaN(e[detJ])) continue; hcoinBGO->Fill(detID, detJ); } - } + + if ( saveEV2){ + + short * out0 = new short[1]; + short * outa = new short[1]; + short * outb = new short[1]; + + out0[0] = numGatedData; + fwrite(out0, 1, 1, outEV2); + + for( int i = 0; i < (int) gatedID.size(); i++){ + int id = gatedID[i]; + outa[0] = id; + fwrite(outa, 1, 1, outEV2); + outb[0] = eCal[id]; + fwrite(outb, 2, 1, outEV2); + } + + fwrite(out0, 1, 1, outEV2); + + /** + int len = (int) gatedID.size(); + char out[2*len+2]; + out[0] = numGatedData; + for( int i = 0; i < (int) gatedID.size(); i++){ + int id = gatedID[i]; + out[2*i+1] = id; + out[2*i+2] = eCal[id]; + } + out[2*len+1]=numGatedData; + fwrite(out,3*len+2,sizeof(out),outEV2); + */ + } + + return kTRUE; } //############################################ TERMINATE void Analyzer::Terminate(){ + + if(saveEV2) fclose(outEV2); printf("============================== finishing.\n"); gROOT->cd(); @@ -177,11 +227,14 @@ void Analyzer::Terminate(){ cCanvas->cd(3); cCanvas->cd(3)->SetLogz(1); - hcoin->Draw("colz"); + hcrystalBGO->Draw("colz"); cCanvas->cd(4); - cCanvas->cd(4)->SetLogz(1); - hcoinBGO->Draw("colz"); + //cCanvas->cd(4)->SetLogz(1); + he[0]->SetLineColor(2); + he[0]->Draw(); + heCal[0]->Draw("same"); + //hcoinBGO->Draw("colz"); printf("=============== loaded AutoFit.C, try showFitMethos()\n"); gROOT->ProcessLine(".L armory/AutoFit.C"); diff --git a/Analyzer.h b/Analyzer.h index b2e900e..1b09294 100644 --- a/Analyzer.h +++ b/Analyzer.h @@ -13,6 +13,7 @@ #include #include #include +#include #include "mapping.h" #include "armory/AnalysisLibrary.h" @@ -28,9 +29,11 @@ public : // Declaration of leaf types ULong64_t evID; Double_t e[NCRYSTAL]; - ULong64_t t[NCRYSTAL]; + ULong64_t e_t[NCRYSTAL]; UShort_t p[NCRYSTAL]; + UShort_t hit[NCRYSTAL]; Double_t bgo[NBGO]; + ULong64_t bgo_t[NBGO]; Double_t other[NOTHER]; Int_t multi; @@ -39,11 +42,18 @@ public : TBranch *b_energy; //! TBranch *b_time; //! TBranch *b_pileup; //! + TBranch *b_hit; //! TBranch *b_bgo; //! + TBranch *b_bgoTime; //! TBranch *b_other; //! - TBranch *b_multiplicity; //! + TBranch *b_multi; //! - Analyzer(TTree * /*tree*/ =0) : fChain(0) { } + Analyzer(TTree * /*tree*/ =0) : fChain(0) { totnumEntry = 0; + Frac = 0.1; + ProcessedEntries = 0; + saveEV2 = true; + outEV2Name = "test.ev2"; + } virtual ~Analyzer() { } virtual Int_t Version() const { return 2; } virtual void Begin(TTree *tree); @@ -60,6 +70,21 @@ public : virtual void Terminate(); ClassDef(Analyzer,0); + + ULong64_t totnumEntry; + + vector> eCorr; + + ULong64_t ProcessedEntries; + Float_t Frac; ///Progress bar + TStopwatch StpWatch; + + bool saveEV2; + FILE * outEV2; + TString outEV2Name; + + + }; #endif @@ -82,11 +107,23 @@ void Analyzer::Init(TTree *tree) fChain->SetBranchAddress("evID", &evID, &b_event_ID); fChain->SetBranchAddress("e", e, &b_energy); - fChain->SetBranchAddress("t", t, &b_time); + fChain->SetBranchAddress("e_t", e_t, &b_time); fChain->SetBranchAddress("p", p, &b_pileup); + fChain->SetBranchAddress("hit", hit, &b_hit); fChain->SetBranchAddress("bgo", bgo, &b_bgo); + fChain->SetBranchAddress("bgo_t", bgo_t, &b_bgoTime); fChain->SetBranchAddress("other", other, &b_other); - fChain->SetBranchAddress("multi", &multi, &b_multiplicity); + fChain->SetBranchAddress("multi", &multi, &b_multi); + + if( saveEV2 ){ + printf("======================== Create output ev2 : %s\n", outEV2Name.Data()); + outEV2 = fopen(outEV2Name.Data(), "w+"); + } + + + printf("======================== Start processing....\n"); + StpWatch.Start(); + } Bool_t Analyzer::Notify() diff --git a/armory/Analyzer_Utili.c b/armory/Analyzer_Utili.c index 5a8f93d..eacc73f 100644 --- a/armory/Analyzer_Utili.c +++ b/armory/Analyzer_Utili.c @@ -18,12 +18,12 @@ void newCanvas(int sizeX = 800, int sizeY = 600, int posX = 0, int posY = 0){ cNewCanvas->cd(); } -void rawEvID(bool cal = false){ - TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID"); - if( cRawID == NULL ) cRawID = new TCanvas("cRawID", "raw ID", 1000, 800); - cRawID->cd(1)->SetGrid(); - cal ? heCalVID->Draw("colz") : heVID->Draw("colz"); -} +//void rawEvID(bool cal = false){ +// TCanvas * cRawID = (TCanvas *) gROOT->FindObjectAny("cRawID"); +// if( cRawID == NULL ) cRawID = new TCanvas("cRawID", "raw ID", 1000, 800); +// cRawID->cd(1)->SetGrid(); +// cal ? heCalVID->Draw("colz") : heVID->Draw("colz"); +//} void drawE(int CloverID = -1, bool cali = false, bool isLogy = false, double xMin = 0, double xMax = 0){ diff --git a/armory/EventBuilder.cpp b/armory/EventBuilder.cpp new file mode 100644 index 0000000..3119d66 --- /dev/null +++ b/armory/EventBuilder.cpp @@ -0,0 +1,203 @@ +#include +#include +#include +#include +#include + +#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] \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 _raw.root \n"); + 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; + + printf("======================== Create output tree\n"); + TString outFileName = inFileName; + outFileName.Remove(inFileName.First('.')); + outFileName.Append(".root"); + if( argn >=4 ) outFileName = argv[3]; + + 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); + +} diff --git a/armory/ev22txt.cpp b/armory/ev22txt.cpp new file mode 100644 index 0000000..7c44366 --- /dev/null +++ b/armory/ev22txt.cpp @@ -0,0 +1,81 @@ +#include +#include +#include +#include + +using namespace std; + +int main(int argc, char **argv){ + + printf("=====================================\n"); + printf("=== ev2 --> txt ===\n"); + printf("=====================================\n"); + + if ( argc != 2 && argc != 3 ){ + printf("Usage: \n"); + printf("%10s [ev2 file] [max number of event]\n", argv[0]); + } + + string inFileName = argv[1]; + string outFileName = argv[1]; + int found = outFileName.find('.'); + int len = outFileName.length(); + outFileName.erase(found, len-found); + outFileName.append(".txt"); + + long long maxNumEvent = 0; // if maxNumEvent <= 0, all event + if( argc >= 3 ) maxNumEvent = atoi(argv[2]); + + + //open list-mode data file from PXI digitizer + FILE *inFile = fopen(argv[1], "r"); + long int inFileSize,inFilepos; + if ( inFile == NULL) { + printf("Error, cannot open input file %s\n", argv[1]); + return 1; + } + + //get file size + fseek(inFile, 0L, SEEK_END); + inFileSize = ftell(inFile); + rewind(inFile); + long int inFilePos = 0; + + printf(" in file : %s | file size : %f MB\n", inFileName.c_str(), inFileSize/1024./1024.); + //printf(" out file : %s \n", outFileName.c_str()); + if( maxNumEvent <= 0 ){ + printf(" max number of event : ALL \n"); + }else{ + printf(" max number of event : %lld \n", maxNumEvent); + } + printf("-------------------------------------\n"); + + long long nEvent = 0; + short * multi = new short[1]; + short * id = new short[1]; + short * energy = new short[1] ; + + while ( inFilePos < inFileSize || feof(inFile) ) { + + printf("------------------- %lld\n", nEvent); + + fread(multi, 1, 1, inFile); + printf(" number of det : %d \n", multi[0]); + + for( int i = 0; i < multi[0] ; i++){ + + fread(id, 1, 1, inFile); + fread(energy, 2, 1, inFile); + + printf("%2d| %3d, %8d \n", i, id[0], energy[0]); + + } + fread(multi, 1, 1, inFile); + + nEvent++; + + if( maxNumEvent > 0 && nEvent > maxNumEvent ) break; + + } + +} diff --git a/armory/evt2root.cpp b/armory/evt2root.cpp index e7b3126..ca5d8b3 100644 --- a/armory/evt2root.cpp +++ b/armory/evt2root.cpp @@ -15,6 +15,7 @@ #define MAX_CHANNELS_PER_BOARD 16 #define BOARD_START 2 +#include "../mapping.h" class measurment{ @@ -106,8 +107,8 @@ int main(int argn, char **argv) { tree->Branch("evID", &measureID, "data_ID/L"); tree->Branch("id", &data.id, "ID/s"); - tree->Branch("e", &data.energy, "energy/s"); - tree->Branch("t", &data.time, "timestamp/l"); + tree->Branch("e", &data.energy, "crystal_energy/s"); + tree->Branch("e_t", &data.time, "crystal_timestamp/l"); //tree->Branch("tdiff", &data.timeDiff, "time_Diff/L"); @@ -125,7 +126,7 @@ int main(int argn, char **argv) { measureID ++; /// see the Pixie-16 user manual, Table4-2 - data.ch = header[0] & 0xF ; + data.ch = header[0] & 0xF ; data.slot = (header[0] >> 4) & 0xF; data.crate = (header[0] >> 8) & 0xF; data.headerLength = (header[0] >> 12) & 0x1F; @@ -138,24 +139,12 @@ int main(int argn, char **argv) { 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.id = mapping[data.id]; + data.energy = data.energy / 2. ; // factor 2 is the rawe_rebin_factor; + nWords += data.eventLength; - //if( measureID == 0 ) { - // data.timeDiff = 0; - //}else{ - // data.timeDiff = (Long64_t) data.time - timeLast; - //} - //timeLast = data.time; - - //if( data.timeDiff == false ){ - // printf("----------------------nWords: %llu, inFilePos: %lu\n", nWords, inFilePos); - // for(int i = 0; i < 4; i++){ - // printf(" %x\n", header[i]); - // } - // data.Print(); - //} - //=== jump to next measurement if( data.eventLength > 4 ){ fseek(inFile, sizeof(int) * (data.eventLength-4), SEEK_CUR); diff --git a/armory/pixie2root.cpp b/armory/pixie2root.cpp deleted file mode 100644 index 1f6c20f..0000000 --- a/armory/pixie2root.cpp +++ /dev/null @@ -1,454 +0,0 @@ -/**********************************************************/ -/* */ -/* Modified by Ryan From */ -/* */ -/* PXI SCAN CODE -- J.M. Allmond (ORNL) -- July 2016 */ -/* */ -/**********************************************************/ - -#include -#include -#include -#include -#include - -#include "TFile.h" -#include "TTree.h" -#include "TMath.h" -#include "TBenchmark.h" - -#define RAND ((float) rand() / ((unsigned int) RAND_MAX + 1)) // random number in interval (0,1) - -#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 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}; - -/**** -int overwrite = 1; - -/////////////////////// -// 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 root ===\n"); - printf("=====================================\n"); - - //CERN ROOT things - TString inFileName = argv[1]; - TString outFileName = inFileName; - - int EVENT_BUILD_TIME = 100; - - if( argc >= 3 ){ - EVENT_BUILD_TIME = atoi(argv[2]); - } - - 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); - printf("------------------------ Event building time window : %d tics = %d nsec \n", EVENT_BUILD_TIME, EVENT_BUILD_TIME*10); - - TFile * outRootFile = new TFile(outFileName, "recreate"); - outRootFile->cd(); - TTree * tree = new TTree("tree", "tree"); - - unsigned long long evID = -1; - double energy[NCRYSTAL]; - unsigned long long etimestamp[NCRYSTAL]; - double bgo[NBGO]; - double other[NOTHER]; - unsigned short pileup[NCRYSTAL]; - - //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", NCRYSTAL)); - tree->Branch("e", energy, Form("energy[%d]/D", NCRYSTAL)); - tree->Branch("t", etimestamp, Form("energy_time_stamp[%d]/l", NCRYSTAL)); - tree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCRYSTAL)); - - 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 < NCRYSTAL; 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 = mapping[subevt[sevtmult].id]; - if ( 0 <= ch && ch < NCRYSTAL ){ - 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 )\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 )\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; -} - - diff --git a/armory/to2root.cpp b/armory/to2root.cpp new file mode 100644 index 0000000..42f01a3 --- /dev/null +++ b/armory/to2root.cpp @@ -0,0 +1,380 @@ +/**********************************************************/ +/* */ +/* Modified by Ryan From */ +/* */ +/* PXI SCAN CODE -- J.M. Allmond (ORNL) -- July 2016 */ +/* */ +/**********************************************************/ + +#include +#include +#include +#include +#include + +#include "TFile.h" +#include "TTree.h" +#include "TMath.h" +#include "TBenchmark.h" + +#define RAND ((float) rand() / ((unsigned int) RAND_MAX + 1)) // random number in interval (0,1) + +#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 RAWE_REBIN_FACTOR 2.0 // Rebin 32k pixie16 spectra to something smaller to fit better into 8k. + +#include "../mapping.h" + +///////////////////// +// RAW EVENT TYPES // +///////////////////// +struct measurement +{ + int chn; + int sln; + int crn; + int id; + int hlen; + int elen; + int trlen; //number of samples + int trwlen; //number of words (two samples per word) + int fcode; //pileup flag + long long int time; + int ctime; + int ctimef; + int e; + int extra; + short tr[4096]; + int esum[4]; + int qsum[8]; +}; +struct measurement data = {0}; + +int sevtmult=0; +unsigned long long int dataCount=0; +unsigned long long int pileUpCount=0; +unsigned long long int evtCount=0; + +/////////////////////////////////// +// START OF MAIN FUNCTION // +/////////////////////////////////// +int main(int argc, char **argv) { + + float tempf=0; + + //temp buffer for each sub event + unsigned int sub[MAX_SUB_LENGTH]; + memset(sub, 0, sizeof(sub)); + + printf("=====================================\n"); + printf("=== evt.to --> root ===\n"); + printf("=====================================\n"); + + // Check that the corrent number of arguments were provided. + if (argc != 2 && argc != 3 ) { + printf("Incorrect number of arguments:\n"); + printf("%s [*.to File] [timeWindow] \n", argv[0]); + printf(" timeWindow : number of tick, 1 tick = 10 ns. default = 100 \n"); + return 1; + } + + //CERN ROOT things + TString inFileName = argv[1]; + TString outFileName = inFileName; + + int timeWindow = 100; + if( argc >= 3 ) timeWindow = atoi(argv[2]); + + outFileName.Remove(inFileName.First('.')); + outFileName.Append(".root"); + + printf(" in file : %s \n", inFileName.Data()); + printf(" our file : %s \n", outFileName.Data()); + + printf(" max number of detector channal: %d \n", MAX_ID); + printf("------------------------ Event building time window : %d tics = %d nsec \n", timeWindow, timeWindow*10); + + TFile * outRootFile = new TFile(outFileName, "recreate"); + outRootFile->cd(); + TTree * tree = new TTree("tree", "tree"); + + unsigned long long evID = -1; + + double e[NCRYSTAL]; + unsigned long long e_t[NCRYSTAL]; + unsigned short pileup[NCRYSTAL]; + unsigned short hit[NCRYSTAL]; // number of hit in an event + + double bgo[NBGO]; + unsigned long long bgo_t[NBGO]; + + double other[NOTHER]; + + int multi; //sum of all crystal hit in an event + + tree->Branch("evID", &evID, "event_ID/l"); + + //TODO: use TCloneArray to save measurement struc, that can save space and possibly time. + ///tree->Branch("detID", detID, Form("det ID[%d]/B", NCRYSTAL)); + tree->Branch("e", e, Form("e[%d]/D", NCRYSTAL)); + tree->Branch("e_t", e_t, Form("e_timestamp[%d]/l", NCRYSTAL)); + tree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCRYSTAL)); + tree->Branch("hit", hit, Form("hit[%d]/s", NCRYSTAL)); + + tree->Branch("bgo", bgo, Form("BGO_e[%d]/D", NBGO)); + tree->Branch("bgo_t", bgo_t, Form("BGO_timestamp[%d]/l", NBGO)); + + tree->Branch("other", other, Form("other_e[%d]/D", NOTHER)); + + tree->Branch("multi", &multi, "multiplicity_crystal/I"); + + //open list-mode data file from PXI digitizer + FILE *fpr = fopen(argv[1], "r"); + long int fprsize,fprpos; + if ( fpr == NULL) { + fprintf(stderr, "Error, cannot open input file %s\n", argv[2]); + return 1; + } + + //get file size + fseek(fpr, 0L, SEEK_END); + fprsize = ftell(fpr); + rewind(fpr); + + TBenchmark gClock; + gClock.Reset(); + gClock.Start("timer"); + + int hitcrystal[NCRYSTAL] = {0}; + + ///////////////////// + // MAIN WHILE LOOP // + ///////////////////// + while (1) { //main while loop + + ///////////////////////////////// + // UNPACK DATA AND EVENT BUILD // + ///////////////////////////////// + + //data clear + 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; + evID++; + + long long int etime=-1; + long long int tdif=-1; + int sevtmult=0; + + while (1) { //get subevents and event build for one "event" + + if (fread(sub, sizeof(int)*HEADER_LENGTH, 1, fpr) != 1) break; + + data.chn = sub[0] & 0xF; /// channel in digitizer + data.sln = (sub[0] & 0xF0) >> 4; /// digitizer ID + data.crn = (sub[0] & 0xF00) >> 8; /// crate + data.id = data.crn*MAX_BOARDS_PER_CRATE*MAX_CHANNELS_PER_BOARD + (data.sln-BOARD_START)*MAX_CHANNELS_PER_BOARD + data.chn; + data.hlen = (sub[0] & 0x1F000) >> 12; + data.elen = (sub[0] & 0x7FFE0000) >> 17; + data.fcode = (sub[0] & 0x80000000) >> 31; + data.time = ( (long long int)(sub[2] & 0xFFFF) << 32) + sub[1]; + data.ctime = (sub[2] & 0x7FFF0000) >> 16; + data.ctimef = (sub[2] & 0x80000000) >> 31; + data.e = (sub[3] & 0xFFFF); + data.trlen = (sub[3] & 0x7FFF0000) >> 16; + data.trwlen = data.trlen / 2; + data.extra = (sub[3] & 0x80000000) >> 31; + + tempf = (float)data.e/RAWE_REBIN_FACTOR;// + RAND; + data.e = (int)tempf; + + //check lengths (sometimes all of the bits for trace length are turned on ...) + /**if (dataBlock[sevtmult].elen - dataBlock[sevtmult].hlen != dataBlock[sevtmult].trwlen) { + printf("SEVERE ERROR: event, header, and trace length inconsistencies found\n"); + printf("event length = %d\n", dataBlock[sevtmult].elen); + printf("header length = %d\n", dataBlock[sevtmult].hlen); + printf("trace length = %d\n", dataBlock[sevtmult].trwlen); + printf("Extra = %d\n", dataBlock[sevtmult].extra); + printf("fcode = %d\n", dataBlock[sevtmult].fcode); + //sleep(1); + //return 0; + } */ + + //Set reference time for event building + if (etime == -1) { + etime = data.time; + tdif = 0; + } + else { + tdif = data.time - etime; + } + + //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! + break; + }else{ + //if within time window, fill array; + int detID = mapping[data.id]; + if ( 0 <= detID && detID < NCRYSTAL ){ + e[detID] = data.e; + e_t[detID] = data.time; + pileup[detID] = data.fcode; + hit[detID] ++; + multi++; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = data.e; + bgo_t[detID-100] = data.time; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = data.e; + } + } + + // total pileups + if (data.fcode==1) { + pileUpCount++; + } + + + //more data than just the header; read entire sub event, first rewind, then read data.elen + fseek(fpr, -sizeof(int)*HEADER_LENGTH, SEEK_CUR); + //if (fread(sub, sizeof(int)*dataBlock[sevtmult].elen, 1, fpr) != 1) break; + if (fread(sub, sizeof(int)*data.elen, 1, fpr) != 1) break; + + /** + //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]; + } + } + */ + 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 + dataCount += sevtmult; + evtCount++; //event-built number + + + int hit_add = 0; + for ( int i = 0; i < NCRYSTAL; i++){ + if( hit[i] > 1 ) { + hit_add = 1; + hitcrystal[i] ++; + } + } + + ///if( hit_add == 1 ){ + /// for ( int i = 0; i < NCRYSTAL; i++){ printf("%d", hit[i]); } + /// printf("------------%d\n", hitcrystal); + ///} + + ///if( evtCount < 40 ) { + /// printf("------------------------------------- %lld \n", evtCount); + /// for( int i = 0; i < NCRYSTAL; i++){ + /// if(e[i] > 0 )printf("%2d %6.0f %10llu\n", i, e[i], e_t[i]); + /// } + /// + ///} + + ///////////////////////////////////// + // END UNPACK DATA AND EVENT BUILD // + ///////////////////////////////////// + + //event stats, print status every 10000 events + if ( evtCount % 10000 == 0 ) { + fprpos = ftell(fpr); + tempf = (float)fprsize/(1024.*1024.*1024.); + gClock.Stop("timer"); + double time = gClock.GetRealTime("timer"); + gClock.Start("timer"); + printf("Total dataBlock: \x1B[32m%llu \x1B[31m(%d%% pileup)\x1B[0m\nTotal Events: \x1B[32m%llu (%.1f )\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\nTime used:%3.0f min %5.2f sec\033[3A\r", + dataCount, (int)((100*pileUpCount)/dataCount), evtCount, (float)dataCount/(float)evtCount, (100*fprpos/fprsize), tempf, TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + } + + + 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 )\x1B[0m\nPercent Complete: \x1B[32m%ld%% of %.3f GB\x1B[0m\n\033[3A\n", + dataCount, (int)((100*pileUpCount)/dataCount), evtCount, (float)dataCount/(float)evtCount, (100*fprpos/fprsize), tempf); + + + outRootFile->cd(); + tree->Write(); + outRootFile->Close(); + + fclose(fpr); + + gClock.Stop("timer"); + double time = gClock.GetRealTime("timer"); + printf("\n\n==================== finished.\r\n"); + printf("Total time spend : %3.0f min %5.2f sec\n", TMath::Floor(time/60.), time - TMath::Floor(time/60.)*60.); + + printf("number of hit per crystal per event:\n"); + for( int i = 0; i < NCRYSTAL ; i++){ + printf("%2d -- %d \n", i, hitcrystal[i]); + } + return 0; +} + + diff --git a/makefile b/makefile index fd541b6..c6537b0 100644 --- a/makefile +++ b/makefile @@ -2,7 +2,7 @@ CC=g++ #all: xia2root xia2ev2_nopart pixie2root scan pxi-time-order #all: xia2root xia2ev2_nopart pixie2root scan evt2root evt2hist -all: xia2root pixie2root evt2root evt2hist pxi-time-order +all: xia2root to2root evt2root evt2hist pxi-time-order ev22txt EventBuilder #this is FSU evt to root xia2root: armory/xia2root.cpp @@ -12,8 +12,8 @@ xia2root: armory/xia2root.cpp # $(CC) armory/xia2ev2_nopart.cpp -o xia2ev2_nopart #this is for eventbuild -pixie2root: armory/pixie2root.cpp - $(CC) armory/pixie2root.cpp -o pixie2root `root-config --cflags --glibs` +to2root: armory/to2root.cpp + $(CC) armory/to2root.cpp -o to2root `root-config --cflags --glibs` #this is for online root evt2root: armory/evt2root.cpp @@ -25,3 +25,9 @@ evt2hist: armory/evt2hist.cpp pxi-time-order: armory/pxi-time-order.c $(CC) armory/pxi-time-order.c -o pxi-time-order + +ev22txt: armory/ev22txt.cpp + $(CC) armory/ev22txt.cpp -o ev22txt + +EventBuilder: armory/EventBuilder.cpp + $(CC) armory/EventBuilder.cpp -o EventBuilder `root-config --cflags --glibs` diff --git a/mapping.h b/mapping.h index d7f4199..b7afe76 100644 --- a/mapping.h +++ b/mapping.h @@ -11,18 +11,18 @@ Other : 200 - 299 #define NBGO 9 #define NOTHER 52 -// 0 1 2 3 4 5 6 7 8 9 -int mapping[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, -1, -1, -1, -1, -1, -1, -1, -1, // 110 - -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // 120 +// 0 1 2 3 4 5 6 7 8 9 +int mapping[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, -1, -1, -1, -1, -1, -1, -1, -1, // 110 + -1, -1, -1, -1, -1, -1, -1, -1, -1, -1}; // 120