Small change of EventBuilder to support TTreeReader. add SpliitPolePlotter.C

This commit is contained in:
Ryan Tang 2024-07-17 18:41:47 -04:00
parent 064c99256b
commit 00ddf3dcf5
5 changed files with 328 additions and 157 deletions

82
.vscode/settings.json vendored
View File

@ -82,86 +82,6 @@
}
],
"files.associations": {
"mainWindow.C": "cpp",
"Scope.C": "cpp",
"new": "cpp",
"allocator": "cpp",
"array": "cpp",
"istream": "cpp",
"ostream": "cpp",
"sstream": "cpp",
"limits": "cpp",
"atomic": "cpp",
"bit": "cpp",
"*.tcc": "cpp",
"bitset": "cpp",
"cctype": "cpp",
"chrono": "cpp",
"clocale": "cpp",
"cmath": "cpp",
"codecvt": "cpp",
"compare": "cpp",
"concepts": "cpp",
"condition_variable": "cpp",
"cstdarg": "cpp",
"cstddef": "cpp",
"cstdint": "cpp",
"cstdio": "cpp",
"cstdlib": "cpp",
"cstring": "cpp",
"ctime": "cpp",
"cwchar": "cpp",
"cwctype": "cpp",
"deque": "cpp",
"list": "cpp",
"map": "cpp",
"string": "cpp",
"unordered_map": "cpp",
"unordered_set": "cpp",
"vector": "cpp",
"exception": "cpp",
"algorithm": "cpp",
"functional": "cpp",
"iterator": "cpp",
"memory": "cpp",
"memory_resource": "cpp",
"numeric": "cpp",
"optional": "cpp",
"random": "cpp",
"ratio": "cpp",
"string_view": "cpp",
"system_error": "cpp",
"tuple": "cpp",
"type_traits": "cpp",
"utility": "cpp",
"future": "cpp",
"initializer_list": "cpp",
"iomanip": "cpp",
"iosfwd": "cpp",
"iostream": "cpp",
"mutex": "cpp",
"numbers": "cpp",
"semaphore": "cpp",
"span": "cpp",
"stdexcept": "cpp",
"stop_token": "cpp",
"streambuf": "cpp",
"thread": "cpp",
"cinttypes": "cpp",
"typeinfo": "cpp",
"variant": "cpp",
"qmainwindow": "cpp",
"qchartview": "cpp",
"qthread": "cpp",
"qrandomgenerator": "cpp",
"source_location": "cpp",
"splitpole.C": "cpp",
"forward_list": "cpp",
"fstream": "cpp",
"Analyzer.C": "cpp",
"process_Run.C": "cpp",
"EncoreAnalyzer.C": "cpp",
"qfiledialog": "cpp",
"script.C": "cpp"
"*.C": "cpp"
}
}

View File

@ -159,8 +159,8 @@ int main(int argc, char **argv) {
tree->Branch("ch", ch, "ch[multi]/s");
tree->Branch("e", e, "e[multi]/s");
tree->Branch("e2", e2, "e2[multi]/s");
tree->Branch("e_t", e_t, "e_timestamp[multi]/l");
tree->Branch("e_f", e_f, "e_fineTime[multi]/s");
tree->Branch("e_t", e_t, "e_t[multi]/l");
tree->Branch("e_f", e_f, "e_f[multi]/s");
tree->Branch("traceLength", traceLength, "traceLength[multi]/s");
if( traceOn ) {

241
Aux/SplitPolePlotter.C Normal file
View File

@ -0,0 +1,241 @@
#include "TFile.h"
#include "TChain.h"
#include "TH1F.h"
#include "TTreeReader.h"
#include "TTreeReaderValue.h"
#include "TTreeReaderArray.h"
#include "TClonesArray.h"
#include "TGraph.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TStopwatch.h"
#include "TMath.h"
#include "vector"
#include "../analyzers/SplitPoleHit.h"
namespace ChMap{
const short ScinR = 0;
const short ScinL = 1;
const short dFR = 8;
const short dFL = 9;
const short dBR = 10;
const short dBL = 11;
const short Cathode = 7;
const short AnodeF = 13;
const short AnodeB = 15;
};
const double c = 299.792458; // mm/ns
const double pi = M_PI;
const double deg2rad = pi/180.;
SplitPoleHit hit;
TH2F * PID;
TH2F * coin;
TH1F * hMulti;
TH1F * hF;
TH1F * hB;
TH1F * hXavg;
TH2F * hFocal;
TH2F * hXavgVQ;
TH2F * haha;
TH1F * hEx;
ULong64_t t1, t2;
#define XMIN -20
#define XMAX 100
//^###########################################
void SplitPolePlotter(TChain *tree){
printf("#####################################################################\n");
printf("################# SplitPolePlotter.C ####################\n");
printf("#####################################################################\n");
TObjArray * fileList = tree->GetListOfFiles();
printf("\033[0;31m========================================== Number of Files : %2d\n",fileList->GetEntries());
fileList->Print();
printf("========================================== Number of Files : %2d\033[0m\n",fileList->GetEntries());
printf("///////////////////////////////////////////////////////////////////\n");
printf(" Total Number of entries : %llu \n", tree->GetEntries());
printf("///////////////////////////////////////////////////////////////////\n");
if( tree->GetEntries() == 0 ) {
printf("========= no events. Abort.\n");
return;
}
//*====================================================== histograms
PID = new TH2F("hPID", "PID; Scin_X ; AnodeB", 200, 0, 30000, 100, 0, 70000);
coin = new TH2F("hCoin", "Coincident ", 16, 0, 16, 16, 0, 16);
hMulti = new TH1F("hMulti", "Multiplicity", 16, 0, 16);
hF = new TH1F("hF", "Front delay line position", 600, XMIN, XMAX);
hB = new TH1F("hB", "Back delay line position", 600, XMIN, XMAX);
hXavg = new TH1F("hAvg", "Xavg", 600, XMIN, XMAX);
hFocal = new TH2F("hFocal", "Front vs Back ", 200, XMIN, XMAX, 200, XMIN, XMAX);
hXavgVQ = new TH2F("hXavgVQ", "Xavg vs Q ", 200, XMIN, XMAX, 200, 0, 40000);
haha = new TH2F("haha", "", 400, XMIN, XMAX, 400, -50, 50);
hEx = new TH1F("hEx", "Ex; Ex [MeV]; count/100 keV", 250, -5, 20);
hit.SetMassTablePath("../analyzers/mass20.txt");
hit.CalConstants("12C", "12C", "4He", 80, 5); // 80MeV, 5 deg
hit.CalZoffset(1.41); // 1.41 T
t1 = 0;
t2 = 0;
//*====================================================== Tree Reader
TTreeReader reader(tree);
TTreeReaderValue<ULong64_t> evID = {reader, "evID"};
TTreeReaderValue<UInt_t> multi = {reader, "multi"};
TTreeReaderArray<UShort_t> sn = {reader, "sn"};
TTreeReaderArray<UShort_t> ch = {reader, "ch"};
TTreeReaderArray<UShort_t> e = {reader, "e"};
TTreeReaderArray<UShort_t> e2 = {reader, "e2"};
TTreeReaderArray<ULong64_t> e_t = {reader, "e_t"};
TTreeReaderArray<UShort_t> e_f = {reader, "e_f"};
ULong64_t NumEntries = tree->GetEntries();
//^###########################################################
//^ * Process
//^###########################################################
printf("############################################### Processing...\n");
fflush(stdout); // flush out any printf
ULong64_t processedEntries = 0;
float Frac = 0.1;
TStopwatch StpWatch;
StpWatch.Start();
while (reader.Next()) {
// if( processedEntries > 10 ) break;
// printf("============== %5llu | multi : %d (%zu) \n", processedEntries, multi.Get()[0], sn.GetSize());
// for( int i = 0; i < multi.Get()[0]; i++ ){
// printf(" %d | %5d %2d %7d %10llu\n", i, sn[i], ch[i], e[i], e_t[i]);
// }
hit.ClearData();
hMulti->Fill(sn.GetSize());
for( int i = 0; i < sn.GetSize(); i++){
t2 = e_t[i];
if( t2 < t1 ) printf("entry %lld-%d, timestamp is not in order. %llu, %llu\n", processedEntries, i, t2, t1);
if( i == 0 ) t1 = e_t[i];
if( ch[i] == ChMap::ScinR ) {hit.eSR = e[i]; hit.tSR = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::ScinL ) {hit.eSL = e[i]; hit.tSL = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::dFR ) {hit.eFR = e[i]; hit.tFR = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::dFL ) {hit.eFL = e[i]; hit.tFL = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::dBR ) {hit.eBR = e[i]; hit.tBR = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::dBL ) {hit.eBL = e[i]; hit.tBL = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::Cathode ) {hit.eCath = e[i]; hit.tCath = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::AnodeF ) {hit.eAF = e[i]; hit.tAF = e_t[i]/4 + e_f[i]/1000.;}
if( ch[i] == ChMap::AnodeB ) {hit.eAB = e[i]; hit.tAB = e_t[i]/4 + e_f[i]/1000.;}
for( int j = i+1; j < sn.GetSize(); j++){
coin->Fill(ch[i], ch[j]);
}
}
unsigned int dQ = hit.eAB; // delta Q
unsigned int Qt = hit.eSR; // total Q
if( Qt > 0 && dQ > 0 ) {
PID->Fill(Qt, dQ);
}
// if( hit.eAF < 50000 ) return kTRUE;
// if( hit.eCath == 0 ) return kTRUE;
// if( hit.eCath > 13000 ) return kTRUE;
hit.CalData();
if( !TMath::IsNaN(hit.x1) || !TMath::IsNaN(hit.x2) ) {
hFocal->Fill(hit.x1, hit.x2);
hF->Fill(hit.x1);
hB->Fill(hit.x2);
hXavg->Fill(hit.xAvg);
hXavgVQ->Fill(hit.xAvg, dQ);
for( int i = 0; i < 400; i++){
double y = -50 + 100/400.*i;
double x = (y/42.8625 + 0.5)* ( hit.x2-hit.x1) + hit.x1;
haha->Fill(x,y);
}
double ex = hit.Rho2Ex( (hit.xAvg/100 + 0.363) );
//if( XMIN < hit.xAvg && hit.xAvg < XMAX) printf("x1 : %6.2f, x2 : %6.2f, xAvg %6.2f cm , ex : %f \n", hit.x1, hit.x2, hit.xAvg, ex);
hEx->Fill(ex);
}
//*============================================ Progress Bar
processedEntries ++;
if (processedEntries >= NumEntries*Frac - 1 ) {
TString msg; msg.Form("%llu", NumEntries/1000);
int len = msg.Sizeof();
printf(" %3.0f%% (%*llu/%llu k) processed in %6.1f sec | expect %6.1f sec\n",
Frac*100, len, processedEntries/1000,NumEntries/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac);
fflush(stdout);
StpWatch.Start(kFALSE);
Frac += 0.1;
}
}
//^###########################################################
//^ * Plot
//^###########################################################
TCanvas * canvas = new TCanvas("cc", "Split-Pole", 1800, 1200);
gStyle->SetOptStat("neiou");
canvas->Divide(3, 3);
canvas->cd(1); PID->Draw("colz");
//canvas->cd(2); coin->Draw("colz");
canvas->cd(2); haha->Draw("colz");
canvas->cd(3); hF->Draw();
canvas->cd(4); hB->Draw();
canvas->cd(5); hXavgVQ->Draw("colz");
canvas->cd(6); hXavg->Draw("colz");
canvas->cd(7); hEx->Draw();
canvas->cd(8); coin->Draw("colz");
canvas->cd(9); canvas->cd(9)->SetLogy(); hMulti->Draw();
}

View File

@ -1,90 +1,101 @@
#include "fsuReader.h"
#include "../MultiBuilder.cpp"
// #include "fsuReader.h"
// #include "../MultiBuilder.cpp"
#include "SplitPolePlotter.C"
void script(){
FSUReader * reader = new FSUReader("~/ExpData/testing/.fsu", 16);
Data * data = reader->GetData();
data->tick2ns = 4;
TChain * chain = new TChain("tree");
reader->ScanNumBlock();
chain->Add("data/temp_002_336_1000.root");
// for( int i = 0; i < 500 ; i++ ) reader->ReadNextBlock(0, 0);
// int ch = 5;
// std::vector<unsigned long long > tList;
// int nEvent = 0;
// for( int i = 0; i < data->TotNumNonPileUpEvents[ch]; i++){
// tList.push_back(data->Timestamp[ch][i]);
// printf("%3d | %d %llu \n", i, data->Energy[ch][i], data->Timestamp[ch][i]);
// nEvent ++;
// }
// std::sort(tList.begin(), tList.end());
// unsigned long long dTime = tList.back() - tList.front();
// double sec = dTime * data->tick2ns / 1e9;
// printf("=========== %llu, %llu = %llu | %f sec | %f Hz\n", tList.back(), tList.front(), dTime, sec, nEvent/sec );
//data->PrintStat(0);
SplitPolePlotter(chain);
data->ClearData();
data->ClearTriggerRate();
//^=====================================================
// FSUReader * reader = new FSUReader("~/ExpData/testing/.fsu", 16);
// Data * data = reader->GetData();
// data->tick2ns = 4;
// reader->ScanNumBlock();
// // for( int i = 0; i < 500 ; i++ ) reader->ReadNextBlock(0, 0);
// // int ch = 5;
// // std::vector<unsigned long long > tList;
// // int nEvent = 0;
// // for( int i = 0; i < data->TotNumNonPileUpEvents[ch]; i++){
// // tList.push_back(data->Timestamp[ch][i]);
// // printf("%3d | %d %llu \n", i, data->Energy[ch][i], data->Timestamp[ch][i]);
// // nEvent ++;
// // }
// // std::sort(tList.begin(), tList.end());
// // unsigned long long dTime = tList.back() - tList.front();
// // double sec = dTime * data->tick2ns / 1e9;
// // printf("=========== %llu, %llu = %llu | %f sec | %f Hz\n", tList.back(), tList.front(), dTime, sec, nEvent/sec );
// //data->PrintStat(0);
MultiBuilder * mb = new MultiBuilder(data, reader->GetDPPType(), 334);
mb->SetTimeWindow(10000);
unsigned long totNumBlock = reader->GetTotNumBlock();
int lastDataIndex = 0;
int lastLoopIndex = 0;
for( unsigned long i = 0; i < 2; i++){
reader->ReadNextBlock();
// int maxDataIndex = 0;
// int maxLoopIndex = 0;
// for( int ch = 0; ch < 16 ; ch++){
// if( data->DataIndex[ch] > maxDataIndex ) maxDataIndex = data->DataIndex[ch];
// if( data->LoopIndex[ch] > maxLoopIndex ) maxLoopIndex = data->LoopIndex[ch];
// }
// if( (maxLoopIndex * MaxNData + maxDataIndex) - (lastLoopIndex * MaxNData + lastDataIndex) > MaxNData * 0.05){
// printf("Agg ID : %lu \n", i );
// data->PrintStat();
// data->PrintAllData();
// mb->BuildEvents();
// mb->PrintAllEvent();
// mb->PrintStat();
// lastDataIndex = maxDataIndex;
// lastLoopIndex = maxLoopIndex;
// }
}
// data->ClearData();
// data->ClearTriggerRate();
data->PrintStat();
data->PrintAllData();
//mb->BuildEvents(true);
// MultiBuilder * mb = new MultiBuilder(data, reader->GetDPPType(), 334);
// mb->SetTimeWindow(10000);
// unsigned long totNumBlock = reader->GetTotNumBlock();
// int lastDataIndex = 0;
// int lastLoopIndex = 0;
mb->BuildEventsBackWard(300);
// for( unsigned long i = 0; i < 2; i++){
mb->PrintAllEvent();
mb->PrintStat();
// reader->ReadNextBlock();
// // int maxDataIndex = 0;
// // int maxLoopIndex = 0;
// // for( int ch = 0; ch < 16 ; ch++){
// // if( data->DataIndex[ch] > maxDataIndex ) maxDataIndex = data->DataIndex[ch];
// // if( data->LoopIndex[ch] > maxLoopIndex ) maxLoopIndex = data->LoopIndex[ch];
// // }
// // if( (maxLoopIndex * MaxNData + maxDataIndex) - (lastLoopIndex * MaxNData + lastDataIndex) > MaxNData * 0.05){
// // printf("Agg ID : %lu \n", i );
// // data->PrintStat();
// // data->PrintAllData();
// // mb->BuildEvents();
// // mb->PrintAllEvent();
// // mb->PrintStat();
// // lastDataIndex = maxDataIndex;
// // lastLoopIndex = maxLoopIndex;
// // }
// }
delete mb;
delete reader;
// data->PrintStat();
// data->PrintAllData();
// //mb->BuildEvents(true);
// mb->BuildEventsBackWard(300);
// mb->PrintAllEvent();
// mb->PrintStat();
// delete mb;
// delete reader;
}

View File

@ -7,7 +7,6 @@
#include <TCanvas.h>
#include <TMath.h>
#include "../analyzers/SplitPoleHit.h"
namespace ChMap{