From d2d7a663961eee5cae7f54d643f96b8563a86fc0 Mon Sep 17 00:00:00 2001 From: "Ryan@WorkStation" Date: Thu, 6 Jan 2022 19:29:26 -0500 Subject: [PATCH] change the EventNuilder, so that tree save dynamic array --- Analyzer.C | 82 ++++++------ Analyzer.h | 33 ++--- armory/Analyzer_Utili.c | 4 +- armory/EventBuilder.cpp | 115 +++++++---------- armory/EventBuilder_seperated.cpp | 206 ++++++++++++++++++++++++++++++ armory/makefile | 2 +- process_run | 8 +- 7 files changed, 309 insertions(+), 141 deletions(-) create mode 100644 armory/EventBuilder_seperated.cpp diff --git a/Analyzer.C b/Analyzer.C index 60d023f..b069217 100644 --- a/Analyzer.C +++ b/Analyzer.C @@ -102,64 +102,62 @@ Bool_t Analyzer::Process(Long64_t entry){ b_energy->GetEntry(entry); b_time->GetEntry(entry); - //b_pileup->GetEntry(entry); - //b_hit->GetEntry(entry); - b_bgo->GetEntry(entry); - b_bgoTime->GetEntry(entry); - b_other->GetEntry(entry); b_multi->GetEntry(entry); + b_multiCry->GetEntry(entry); + b_detID->GetEntry(entry); if( multi == 0 ) return kTRUE; for( int i = 0; i < NCRYSTAL; i++) eCal[i] = TMath::QuietNaN(); - ///=========== Looping Crystals - for( int detID = 0; detID < NCRYSTAL ; detID ++){ - - //======== baics gate when no energy or pileup - if( TMath::IsNaN(e[detID])) continue; - //if( pileup[detID] == 1 ) continue; - + ///=========== Looping data for the event + for( int i = 0; i < multi ; i ++){ + //======== Fill raw data - 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; - //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 && 4*kk <= detID && detID < 4*(kk+1) ) { - return kTRUE; + if( detID[i] < 100 ){ /// gamma data + heVID->Fill( detID[i], e[i]); + he[detID[i]]->Fill(e[i]); + + for ( int j = i + 1; j < multi; j++){ + if( 100 <= detID[j] && detID[j] < 200 ) hcrystalBGO->Fill(detID[i], detID[j]-100); /// crystal - BGO coincident + + if( detID[j] < 100 ) hcoin->Fill(detID[i], detID[j]); /// crystal-crystal coincident + } } + if ( 100 < detID[i] && detID[i] < 200 ){ /// BGO data + + } + + + //======== BGO veto + bool dropflag = false; + if( detID[i] < 100 && multi > 1) { + for( int j = i + 1; j < multi; j++){ + if( detID[j] > 100 && (detID[j]-100)*4 < detID[i] && detID[i] < (detID[j]-100 +1)*4) { + dropflag = true; + break; + } + } + } + if( dropflag ) return kTRUE; + //========= 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); //} - if( e_corr == "" ){ - eCal[detID] = e[detID]; - }else{ - eCal[detID] = eCorr[detID][0] + eCorr[detID][1] * e[detID]; - } - 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( detID[i] < 100 ) { + if( e_corr == "" ){ + eCal[detID[i]] = e[i]; + }else{ + eCal[detID[i]] = eCorr[detID[i]][0] + eCorr[detID[i]][1] * e[i]; + } + + heCalVID->Fill( detID[i], eCal[detID[i]]); + heCal[detID[i]]->Fill(eCal[detID[i]]); } } diff --git a/Analyzer.h b/Analyzer.h index 50a9b83..63f17e1 100644 --- a/Analyzer.h +++ b/Analyzer.h @@ -28,25 +28,19 @@ public : // Declaration of leaf types ULong64_t evID; - Double_t e[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 detID[200]; + Double_t e[200]; + ULong64_t e_t[200]; Int_t multi; + Int_t multiCry; // List of branches TBranch *b_event_ID; //! + TBranch *b_detID; //! TBranch *b_energy; //! TBranch *b_time; //! - TBranch *b_pileup; //! - TBranch *b_hit; //! - TBranch *b_bgo; //! - TBranch *b_bgoTime; //! - TBranch *b_other; //! TBranch *b_multi; //! + TBranch *b_multiCry; //! Analyzer(TTree * /*tree*/ =0) : fChain(0) { totnumEntry = 0; Frac = 0.1; @@ -104,15 +98,12 @@ void Analyzer::Init(TTree *tree) fChain = tree; fChain->SetMakeClass(1); - fChain->SetBranchAddress("evID", &evID, &b_event_ID); - fChain->SetBranchAddress("e", e, &b_energy); - 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_multi); + fChain->SetBranchAddress("evID", &evID, &b_event_ID); + fChain->SetBranchAddress("id", detID, &b_detID); + fChain->SetBranchAddress("e", e, &b_energy); + fChain->SetBranchAddress("e_t", e_t, &b_time); + fChain->SetBranchAddress("multi", &multi, &b_multi); + fChain->SetBranchAddress("multiCry", &multiCry, &b_multiCry); TString option = GetOption(); if ( option != "" ) outEV2Name = option; diff --git a/armory/Analyzer_Utili.c b/armory/Analyzer_Utili.c index 99a9cdd..cc3f21c 100644 --- a/armory/Analyzer_Utili.c +++ b/armory/Analyzer_Utili.c @@ -70,10 +70,10 @@ void drawE(int CloverID = -1, bool cali = false, bool isLogy = false, double xMi int hID = nCrystalPerClover * CloverID + i ; if( cali ) { heCal[hID]->SetMaximum(maxY); - heCal[hID]->Draw(""); + heCal[hID]->Draw("same"); }else{ he[hID]->SetMaximum(maxY); - he[hID]->Draw(""); + he[hID]->Draw("same"); } diff --git a/armory/EventBuilder.cpp b/armory/EventBuilder.cpp index 45fd3cf..96add19 100644 --- a/armory/EventBuilder.cpp +++ b/armory/EventBuilder.cpp @@ -13,31 +13,7 @@ #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; -} +#define MAXMULTI 200 int main(int argn, char **argv){ printf("=====================================\n"); @@ -105,26 +81,26 @@ int main(int argn, char **argv){ saveFile->cd(); TTree * newtree = new TTree("tree", "tree"); + Int_t multi = 0; /// this is total multipicilty for all detectors + newtree->Branch("multi", &multi, "multipiclity/I"); + + Int_t eventID = 0 ; 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(); + Int_t multiCry = 0 ; /// thi is total multiplicity for all crystal + newtree->Branch("multiCry", &multiCry, "multiplicity_crystal/I"); + + int id[MAXMULTI]; + double e[MAXMULTI]; + ULong64_t e_t[MAXMULTI]; + newtree->Branch("id", id, "id[multipiclity]/I" ); + newtree->Branch("e", e, "e[multipiclity]/D" ); + newtree->Branch("e_t", e_t, "e_timestamp[multipiclity]/l"); 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++){ @@ -145,52 +121,47 @@ int main(int argn, char **argv){ b_energy->GetEntry(entry); b_energy_timestamp->GetEntry(entry); - if( time0 == 0) time0 = energy_t; + if( time0 == 0) { + time0 = energy_t; + multi = 0; + } 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); + id[multi] = detID; + e[multi] = energy; + e_t[multi] = energy_t; + multi ++; + if( detID < NCRYSTAL ) multiCry++; }else{ - //---- end of event - eventID ++; - + ///---- end of event saveFile->cd(); newtree->Fill(); + eventID ++; + + ///---- clear data + for( int i = 0; i < multi ; i ++){ + id[i] = 0; + e[i] = TMath::QuietNaN(); + e_t[i] = 0; + } + multi = 0; + multiCry = 0; - 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; - } - + + id[multi] = detID; + e[multi] = energy; + e_t[multi] = energy_t; + multi++; + if( detID < NCRYSTAL ) multiCry++; + } } diff --git a/armory/EventBuilder_seperated.cpp b/armory/EventBuilder_seperated.cpp new file mode 100644 index 0000000..45fd3cf --- /dev/null +++ b/armory/EventBuilder_seperated.cpp @@ -0,0 +1,206 @@ +#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 %s \n", inFileName.Data()); + TFile * inFile = new TFile(inFileName, "READ"); + if( inFile->IsOpen() == false ) { + printf("!!!! cannot open file %s \n", inFileName.Data()); + return 0; + } + + TTree * tree = (TTree *) inFile->Get("tree"); + + Long64_t evID; + UShort_t detID; + UShort_t energy; + ULong64_t energy_t; + + TBranch *b_data_ID; //! + TBranch *b_ID; //! + TBranch *b_energy; //! + TBranch *b_energy_timestamp; //! + + tree->SetBranchAddress("evID", &evID, &b_data_ID); + tree->SetBranchAddress("id", &detID, &b_ID); + tree->SetBranchAddress("e", &energy, &b_energy); + tree->SetBranchAddress("e_t", &energy_t, &b_energy_timestamp); + + Long64_t totnumEntry = tree->GetEntries(); + + printf( "total Entry : %lld \n", totnumEntry); + + printf(">>> Buidling Index using the timestamp\n"); + tree->BuildIndex("e_t"); + TTreeIndex *in = (TTreeIndex*) tree->GetTreeIndex(); + Long64_t * index = in->GetIndex(); + + ULong64_t time0; //time-0 for each event + int timeDiff; + + + TString outFileName = inFileName; + outFileName.Remove(inFileName.First("_raw")); + outFileName.Append(".root"); + if( argn >=4 ) outFileName = argv[3]; + + printf(">>> out File name : %s\n", outFileName.Data()); + + printf(">>> Create output tree\n"); + TFile * saveFile = new TFile(outFileName, "recreate"); + saveFile->cd(); + TTree * newtree = new TTree("tree", "tree"); + + newtree->Branch("evID", &eventID, "event_ID/l"); + newtree->Branch("e", e, Form("e[%d]/D", NCRYSTAL)); + newtree->Branch("e_t", e_t, Form("e_timestamp[%d]/l", NCRYSTAL)); + //newtree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCRYSTAL)); + //newtree->Branch("hit", hit, Form("hit[%d]/s", NCRYSTAL)); + + newtree->Branch("bgo", bgo, Form("BGO_e[%d]/D", NBGO)); + newtree->Branch("bgo_t", bgo_t, Form("BGO_timestamp[%d]/l", NBGO)); + + newtree->Branch("other", other, Form("other_e[%d]/D", NOTHER)); + + newtree->Branch("multi", &multi, "multiplicity_crystal/I"); + + ClearTreeData(); + + printf("================== Start processing....\n"); + Float_t Frac = 0.1; ///Progress bar + TStopwatch StpWatch; + StpWatch.Start(); + eventID = 0; + + for( Long64_t entry = 0; entry < totnumEntry; entry++){ + + + /*********** Progress Bar ******************************************/ + if (entry>totnumEntry*Frac-1) { + TString msg; msg.Form("%llu", totnumEntry/1000); + int len = msg.Sizeof(); + printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n", + Frac*100, len, entry/1000,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); + StpWatch.Start(kFALSE); + Frac+=0.1; + } + + entry = index[entry]; + + b_ID->GetEntry(entry); + b_energy->GetEntry(entry); + b_energy_timestamp->GetEntry(entry); + + if( time0 == 0) time0 = energy_t; + timeDiff = (int) (energy_t - time0); + + if( timeDiff < timeWindow ) { + + if ( detID < NCRYSTAL ){ + e[detID] = energy; + e_t[detID] = energy_t; + multi++; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = energy; + bgo_t[detID-100] = energy_t; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = energy; + } + + //printf("%d | %3d %6d %10llu, %3d\n", multi, detID, energy, energy_t, timeDiff); + + }else{ + //---- end of event + eventID ++; + + saveFile->cd(); + newtree->Fill(); + + ClearTreeData(); + + /// fill 1st data of an event + time0 = energy_t; + timeDiff = 0; + + if ( detID < NCRYSTAL ){ + e[detID] = energy; + e_t[detID] = energy_t; + multi = 1; + } + if ( 100 <= detID && detID < 100 + NBGO ){ + bgo[detID-100] = energy; + bgo_t[detID-100] = energy_t; + } + if ( 200 <= detID && detID < 200 + NOTHER){ + other[detID-200] = energy; + } + + } + + } + + printf("============================== finished.\n"); + + saveFile->cd(); + newtree->Write(); + saveFile->Close(); + + printf(" total number of event Built : %d \n", eventID); + +} diff --git a/armory/makefile b/armory/makefile index bba2edb..c66592e 100644 --- a/armory/makefile +++ b/armory/makefile @@ -4,7 +4,7 @@ CC=g++ #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 MergeEVT ev22txt EventBuilder pxi-time-order EventBuilder2 #this is FSU evt to root xia2root: ../armory/xia2root.cpp diff --git a/process_run b/process_run index 94fff15..c21d138 100755 --- a/process_run +++ b/process_run @@ -3,6 +3,8 @@ DIR=$(pwd) DATA_DIR=data +TIMEWINDOW=100 + if [ $# -eq 0 ] || [ $1 == "-help" ]; then echo "$./process_run [Run Folder] [Merge] [isBuildEvents] [Analysis]" echo " Run Folder = the name of run folder" @@ -77,12 +79,12 @@ if [ ${isBuildEvents} -eq 1 ]; then rootDateTime=`stat -c "%Z" ${RunFolder}".root" | sort -rn | head -1` if [ ${rawRootDateTime} -gt ${rootDateTime} ]; then - ./armory/EventBuilder ${RunFolder}"_raw.root" + ./armory/EventBuilder ${RunFolder}"_raw.root" $TIMEWINDOW else echo -e "${RunFolder}.root is up-to-date." fi else - ./armory/EventBuilder ${RunFolder}"_raw.root" + ./armory/EventBuilder ${RunFolder}"_raw.root" $TIMEWINDOW fi fi @@ -93,7 +95,7 @@ fi if [ ${isBuildEvents} -eq -1 ]; then echo -e "$YELLOW forced by user $NC" - ./armory/EventBuilder ${RunFolder}"_raw.root" + ./armory/EventBuilder ${RunFolder}"_raw.root" $TIMEWINDOW fi echo -e "$RED>>> `date` >>>>>>>>>>>>>>>>>>>>>>> Build Events finsihed.$NC"