change the EventNuilder, so that tree save dynamic array

This commit is contained in:
Ryan Tang 2022-01-06 19:29:26 -05:00
parent 273342954a
commit d2d7a66396
7 changed files with 309 additions and 141 deletions

View File

@ -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]);
if( detID[i] < 100 ){ /// gamma data
heVID->Fill( detID[i], e[i]);
he[detID[i]]->Fill(e[i]);
for( int kk = 0 ; kk < NBGO; kk++){
if( bgo[kk] > 0 ) hcrystalBGO->Fill(detID, kk);
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
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;
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( detID[i] < 100 ) {
if( e_corr == "" ){
eCal[detID] = e[detID];
eCal[detID[i]] = e[i];
}else{
eCal[detID] = eCorr[detID][0] + eCorr[detID][1] * e[detID];
eCal[detID[i]] = eCorr[detID[i]][0] + eCorr[detID[i]][1] * e[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);
heCalVID->Fill( detID[i], eCal[detID[i]]);
heCal[detID[i]]->Fill(eCal[detID[i]]);
}
}

View File

@ -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;
@ -105,14 +99,11 @@ void Analyzer::Init(TTree *tree)
fChain->SetMakeClass(1);
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("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("multiCry", &multiCry, &b_multiCry);
TString option = GetOption();
if ( option != "" ) outEV2Name = option;

View File

@ -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");
}

View File

@ -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));
Int_t multiCry = 0 ; /// thi is total multiplicity for all crystal
newtree->Branch("multiCry", &multiCry, "multiplicity_crystal/I");
newtree->Branch("other", other, Form("other_e[%d]/D", NOTHER));
newtree->Branch("multi", &multi, "multiplicity_crystal/I");
ClearTreeData();
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,51 +121,46 @@ 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;
id[multi] = detID;
e[multi] = energy;
e_t[multi] = 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);
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++;
}

View File

@ -0,0 +1,206 @@
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <stdbool.h>
#include "TFile.h"
#include "TTree.h"
#include "TMath.h"
#include "TBenchmark.h"
#include "TStopwatch.h"
#include "TTreeIndex.h"
#include "../mapping.h"
Int_t eventID = 0 ;
double e[NCRYSTAL];
ULong64_t e_t[NCRYSTAL];
double bgo[NBGO];
ULong64_t bgo_t[NBGO];
Short_t other[NOTHER];
Short_t multi;
void ClearTreeData(){
for( int i = 0; i < NCRYSTAL; i++){
e[i] = TMath::QuietNaN();
e_t[i] = 0;
//pileup[i] = 0;
//hit[i] = 0;
}
for( int i = 0; i < NBGO; i++) {
bgo[i] = TMath::QuietNaN();
bgo_t[i] = 0 ;
}
for( int i = 0; i < NOTHER; i++) {
other[i] = TMath::QuietNaN();
}
multi = 0;
}
int main(int argn, char **argv){
printf("=====================================\n");
printf("=== Event Builder ===\n");
printf("=====================================\n");
if (argn != 2 && argn != 3 && argn != 4 ) {
printf("Usage :\n");
printf("%s [_raw.root File] <timeWindows> <SaveFileName>\n", argv[0]);
printf(" timeWindows : default = 100 \n");
printf(" SaveFileName : default is *.root \n");
return 1;
}
TString inFileName = argv[1]; // need to check name
int timeWindow = 100;
if( argn >= 3 ) timeWindow = atoi(argv[2]);
printf(">>> Opening input %s \n", inFileName.Data());
TFile * inFile = new TFile(inFileName, "READ");
if( inFile->IsOpen() == false ) {
printf("!!!! cannot open file %s \n", inFileName.Data());
return 0;
}
TTree * tree = (TTree *) inFile->Get("tree");
Long64_t evID;
UShort_t detID;
UShort_t energy;
ULong64_t energy_t;
TBranch *b_data_ID; //!
TBranch *b_ID; //!
TBranch *b_energy; //!
TBranch *b_energy_timestamp; //!
tree->SetBranchAddress("evID", &evID, &b_data_ID);
tree->SetBranchAddress("id", &detID, &b_ID);
tree->SetBranchAddress("e", &energy, &b_energy);
tree->SetBranchAddress("e_t", &energy_t, &b_energy_timestamp);
Long64_t totnumEntry = tree->GetEntries();
printf( "total Entry : %lld \n", totnumEntry);
printf(">>> Buidling Index using the timestamp\n");
tree->BuildIndex("e_t");
TTreeIndex *in = (TTreeIndex*) tree->GetTreeIndex();
Long64_t * index = in->GetIndex();
ULong64_t time0; //time-0 for each event
int timeDiff;
TString outFileName = inFileName;
outFileName.Remove(inFileName.First("_raw"));
outFileName.Append(".root");
if( argn >=4 ) outFileName = argv[3];
printf(">>> out File name : %s\n", outFileName.Data());
printf(">>> Create output tree\n");
TFile * saveFile = new TFile(outFileName, "recreate");
saveFile->cd();
TTree * newtree = new TTree("tree", "tree");
newtree->Branch("evID", &eventID, "event_ID/l");
newtree->Branch("e", e, Form("e[%d]/D", NCRYSTAL));
newtree->Branch("e_t", e_t, Form("e_timestamp[%d]/l", NCRYSTAL));
//newtree->Branch("p", pileup, Form("pile_up_flag[%d]/s", NCRYSTAL));
//newtree->Branch("hit", hit, Form("hit[%d]/s", NCRYSTAL));
newtree->Branch("bgo", bgo, Form("BGO_e[%d]/D", NBGO));
newtree->Branch("bgo_t", bgo_t, Form("BGO_timestamp[%d]/l", NBGO));
newtree->Branch("other", other, Form("other_e[%d]/D", NOTHER));
newtree->Branch("multi", &multi, "multiplicity_crystal/I");
ClearTreeData();
printf("================== Start processing....\n");
Float_t Frac = 0.1; ///Progress bar
TStopwatch StpWatch;
StpWatch.Start();
eventID = 0;
for( Long64_t entry = 0; entry < totnumEntry; entry++){
/*********** Progress Bar ******************************************/
if (entry>totnumEntry*Frac-1) {
TString msg; msg.Form("%llu", totnumEntry/1000);
int len = msg.Sizeof();
printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n",
Frac*100, len, entry/1000,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac);
StpWatch.Start(kFALSE);
Frac+=0.1;
}
entry = index[entry];
b_ID->GetEntry(entry);
b_energy->GetEntry(entry);
b_energy_timestamp->GetEntry(entry);
if( time0 == 0) time0 = energy_t;
timeDiff = (int) (energy_t - time0);
if( timeDiff < timeWindow ) {
if ( detID < NCRYSTAL ){
e[detID] = energy;
e_t[detID] = energy_t;
multi++;
}
if ( 100 <= detID && detID < 100 + NBGO ){
bgo[detID-100] = energy;
bgo_t[detID-100] = energy_t;
}
if ( 200 <= detID && detID < 200 + NOTHER){
other[detID-200] = energy;
}
//printf("%d | %3d %6d %10llu, %3d\n", multi, detID, energy, energy_t, timeDiff);
}else{
//---- end of event
eventID ++;
saveFile->cd();
newtree->Fill();
ClearTreeData();
/// fill 1st data of an event
time0 = energy_t;
timeDiff = 0;
if ( detID < NCRYSTAL ){
e[detID] = energy;
e_t[detID] = energy_t;
multi = 1;
}
if ( 100 <= detID && detID < 100 + NBGO ){
bgo[detID-100] = energy;
bgo_t[detID-100] = energy_t;
}
if ( 200 <= detID && detID < 200 + NOTHER){
other[detID-200] = energy;
}
}
}
printf("============================== finished.\n");
saveFile->cd();
newtree->Write();
saveFile->Close();
printf(" total number of event Built : %d \n", eventID);
}

View File

@ -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

View File

@ -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"