From 9700accb0f30ac395f505a23ecf1486865712472 Mon Sep 17 00:00:00 2001 From: "nrb@Debain10" Date: Fri, 18 Mar 2022 18:34:13 -0400 Subject: [PATCH] added GAGGPIDCutCreator, PreAnalyzer to reduce the root size, PIDAnalyzer --- GAGGPIDCutCreator.C | 97 ++++++++++++++++ PIDAnalyzer.C | 263 ++++++++++++++++++++++++++++++++++++++++++++ PIDAnalyzer.h | 138 +++++++++++++++++++++++ PreAnalyzer.C | 195 ++++++++++++++++++++++++++++++++ PreAnalyzer.h | 200 +++++++++++++++++++++++++++++++++ 5 files changed, 893 insertions(+) create mode 100644 GAGGPIDCutCreator.C create mode 100644 PIDAnalyzer.C create mode 100644 PIDAnalyzer.h create mode 100644 PreAnalyzer.C create mode 100644 PreAnalyzer.h diff --git a/GAGGPIDCutCreator.C b/GAGGPIDCutCreator.C new file mode 100644 index 0000000..a62837e --- /dev/null +++ b/GAGGPIDCutCreator.C @@ -0,0 +1,97 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "mapping.h" + +void GAGGPIDCutCreator(TString dataList, + TString saveFileName = "GAGGPIDCuts.root", + int peakRange=1200, + int tailRange=400, + bool isLogz = false, + TString gate = "", + TString treeName = "tree" + ){ + + printf("================ Graphic Cut Creator for GAGG ============== \n"); + + TChain * chain = new TChain(treeName); + chain->Add(dataList); + + chain->GetListOfFiles()->Print(); + + TString varX, varY, tag; + + gStyle->SetOptStat("neiou"); + + TCanvas * cCutCreator = new TCanvas("cCutCreator", "GAGG PID Cut Creator", 100, 100, 800, 800); + if( !cCutCreator->GetShowToolBar() ) cCutCreator->ToggleToolBar(); + + cCutCreator->Update(); + if( isLogz ) cCutCreator->cd()->SetLogz(); + + TCutG * cut = NULL; + TObjArray * cutList = new TObjArray(); + + TString expression; + + TH2F * h[NGAGG]; + + int count = 0; + for (Int_t i = 0; i < NGAGG; i++) { + + + varX.Form("gaggT"); varY.Form("gaggP"); + + h[i] = new TH2F(Form("h%d", i), Form("GAGG-%d",i), 500, 0, tailRange, 500, 0, peakRange); + + expression.Form("%s:%s>>h%d", varY.Data(), varX.Data(),i); + gate.Form("gaggID==%d", i); + + chain->Draw(expression, gate, "col"); + + if( h[i]->Integral() < 1000 ) { + h[i]->SetMarkerStyle(20); + h[i]->SetMarkerSize(0.4); + h[i]->Draw(""); + } + + printf("======== make a graphic cut on the plot (double click to stop), %d-th cut: ", i ); + + cCutCreator->Modified(); cCutCreator->Update(); + + gPad->WaitPrimitive(); + + cut = (TCutG*) gROOT->FindObject("CUTG"); + + if( cut == NULL ){ + printf(" stopped by user. no file saved or changed. \n"); + break; + } + + TString name; name.Form("cut%d", i); + cut->SetName(name); + cut->SetVarX(varX.Data()); + cut->SetVarY(varY.Data()); + cut->SetTitle(tag); + cut->SetLineColor(i+1); + cutList->Add(cut); + + printf(" cut-%d \n", i); + + count ++; + } + + TFile * cutFile = new TFile(saveFileName, "recreate"); + cutList->Write("cutList", TObject::kSingleKey); + + printf("====> saved %d cuts into %s\n", count, saveFileName.Data()); + gROOT->ProcessLine(".q"); + +} diff --git a/PIDAnalyzer.C b/PIDAnalyzer.C new file mode 100644 index 0000000..3c0eeee --- /dev/null +++ b/PIDAnalyzer.C @@ -0,0 +1,263 @@ +#define PIDAnalyzer_cxx + + +#include "PIDAnalyzer.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include + +//############################################ User setting + +int energyRange[3] = {1, 30, 800}; // keV {resol, min, max} + +int pidMaxRange[3] = {500, 400, 1600}; //nBin, tail, and peak + +TString cutFileName1 = "testCuts.root"; +TString cutFileName2 = ""; //"alphaCut.root"; +TString cutFileName3 = ""; //"tritonCut.root"; + +short timeGateFlag = 4; // 0 = off, 1 <, 2 >, 3 sandwish, 4 !sandwish +unsigned int timeGate[2] = {45, 65}; // if timeGateFlag < 3, only timeGate[0] use, else, {min, max} + +//############################################ end of user setting + +TH2F * hgg; + +TH1F * hg[NCLOVER]; +TH1F * hg_g1[NCLOVER]; +TH1F * hg_g2[NCLOVER]; +TH1F * hg_g3[NCLOVER]; + +TH2F * hPID[NGAGG]; + +///============= cut +TObjArray * cutList1; +TObjArray * cutList2; +TObjArray * cutList3; +TCutG * cut; + +void PIDAnalyzer::Begin(TTree *tree){ + + TString option = GetOption(); + + totnumEntry = tree->GetEntries(); + + printf( "=========================================================================== \n"); + printf( "========================== Analysis.C/h ================================ \n"); + printf( "====== total Entry : %lld \n", totnumEntry); + printf( "=========================================================================== \n"); + + printf("======================== Histograms declaration\n"); + + + for( int i = 0; i < NCLOVER; i++){ + hg[i] = new TH1F(Form("hg%02d", i), Form("Clover-%02d (added-back)", i), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + hg_g1[i] = new TH1F(Form("hg_g1%02d", i), Form("Clover-%02d (added-back) particle", i), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + hg_g2[i] = new TH1F(Form("hg_g2%02d", i), Form("Clover-%02d (added-back) particle", i), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + hg_g3[i] = new TH1F(Form("hg_g3%02d", i), Form("Clover-%02d (added-back) particle", i), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + } + + hgg = new TH2F("hgg", "Gamma - Gamma", (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2], (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]); + + + for( int i = 0; i < NGAGG; i++){ + hPID[i] = new TH2F(Form("hPID%02d", i), Form("PID-%2d; tail; peak ", i), pidMaxRange[0], -20, pidMaxRange[1], pidMaxRange[0], -50, pidMaxRange[2]); + } + + printf("======================== Load cuts.\n"); + if( cutFileName1 != ""){ + TFile * cutFile1 = new TFile(cutFileName1); + if( cutFile1->IsOpen()){ + cutList1 = (TObjArray *) cutFile1->FindObjectAny("cutList"); + int numCut1 = cutList1->GetEntries(); + printf("=========== found %d cutG in %s \n", numCut1, cutFile1->GetName()); + + for(int i = 0; i < numCut1 ; i++){ + printf("cut name : %s , VarX: %s, VarY: %s, numPoints: %d \n", + cutList1->At(i)->GetName(), + ((TCutG*)cutList1->At(i))->GetVarX(), + ((TCutG*)cutList1->At(i))->GetVarY(), + ((TCutG*)cutList1->At(i))->GetN()); + } + }else{ + cutList1 = NULL; + } + } + + if( cutFileName2 != ""){ + TFile * cutFile2 = new TFile(cutFileName2); + if( cutFile2->IsOpen()){ + cutList2 = (TObjArray *) cutFile2->FindObjectAny("cutList"); + int numCut2 = cutList2->GetEntries(); + printf("=========== found %d cutG in %s \n", numCut2, cutFile2->GetName()); + + for(int i = 0; 2 < numCut2 ; i++){ + printf("cut name : %s , VarX: %s, VarY: %s, numPoints: %d \n", + cutList2->At(i)->GetName(), + ((TCutG*)cutList2->At(i))->GetVarX(), + ((TCutG*)cutList2->At(i))->GetVarY(), + ((TCutG*)cutList2->At(i))->GetN()); + } + }else{ + cutList2 = NULL; + } + } + + if( cutFileName3 != ""){ + TFile * cutFile3 = new TFile(cutFileName3); + if( cutFile3->IsOpen()){ + cutList3 = (TObjArray *) cutFile3->FindObjectAny("cutList"); + int numCut3 = cutList3->GetEntries(); + printf("=========== found %d cutG in %s \n", numCut3, cutFile3->GetName()); + + for(int i = 0; 2 < numCut3 ; i++){ + printf("cut name : %s , VarX: %s, VarY: %s, numPoints: %d \n", + cutList3->At(i)->GetName(), + ((TCutG*)cutList3->At(i))->GetVarX(), + ((TCutG*)cutList3->At(i))->GetVarY(), + ((TCutG*)cutList3->At(i))->GetN()); + } + }else{ + cutList3 = NULL; + } + } + +} + +Bool_t PIDAnalyzer::Process(Long64_t entry){ + ProcessedEntries++; + + /*********** Progress Bar ******************************************/ + 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,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); + StpWatch.Start(kFALSE); + Frac+=0.1; + } + + b_eventID->GetEntry(entry); + b_runID->GetEntry(entry); + b_multi->GetEntry(entry); + b_multiGagg->GetEntry(entry); + b_gammaID->GetEntry(entry); + b_gamma->GetEntry(entry); + b_gamma_t->GetEntry(entry); + b_gaggID->GetEntry(entry); + b_gaggP->GetEntry(entry); + b_gaggT->GetEntry(entry); + b_gagg_t->GetEntry(entry); + + ///################## Gagg + bool fillFlag1 = false; + bool fillFlag2 = false; + bool fillFlag3 = false; + + for( int i = 0 ; i < multiGagg; i++){ + hPID[gaggID[i]]->Fill(gaggT[i],gaggP[i]); + + if( cutList1 != NULL && i < cutList1->GetEntries()) { + cut = (TCutG *) cutList1->At(i); + if( fillFlag1 == false && cut->IsInside(gaggT[i],gaggP[i]) ) { + fillFlag1 = true; + break; + } + } + if( cutList2 != NULL && i < cutList2->GetEntries()) { + cut = (TCutG *) cutList2->At(i); + if( fillFlag2 == false && cut->IsInside(gaggT[i],gaggP[i]) ) { + fillFlag2 = true; + break; + } + } + if( cutList3 != NULL && i < cutList3->GetEntries()) { + cut = (TCutG *) cutList3->At(i); + if( fillFlag2 == false && cut->IsInside(gaggT[i],gaggP[i]) ) { + fillFlag3 = true; + break; + } + } + + } + + + ///################## Gamma data from Clover + + for( int i = 0; i < multi ; i ++){ + + //======== TDiff veto + //if( timeGateFlag == 1 && e_t[i] - e_t[0] > timeGate[0] ) continue; + //if( timeGateFlag == 2 && e_t[i] - e_t[0] < timeGate[0] ) continue; + //if( timeGateFlag == 3 && !(timeGate[0] < e_t[i] - e_t[0] && e_t[i] - e_t[0] < timeGate[1]) ) continue; + //if( timeGateFlag == 4 && timeGate[0] < e_t[i] - e_t[0] && e_t[i] - e_t[0] < timeGate[1] ) continue; + + hg[gammaID[i]]->Fill(gamma[i]); + + if( fillFlag1 ) hg_g1[gammaID[i]]->Fill(gamma[i]); + if( fillFlag2 ) hg_g2[gammaID[i]]->Fill(gamma[i]); + if( fillFlag3 ) hg_g3[gammaID[i]]->Fill(gamma[i]); + + for( int j = 0 ; j < multi; j++){ + if( j != i ) hgg->Fill( gamma[i], gamma[j]); + } + + } + + + return kTRUE; +} + + + +void PIDAnalyzer::Terminate(){ + + printf("============================== finishing.\n"); + gROOT->cd(); + + int canvasXY[2] = {1600 , 800} ;// x, y + int canvasDiv[2] = {4,3}; + TCanvas *cCanvas = new TCanvas("cCanvas", "" ,canvasXY[0],canvasXY[1]); + if( !cCanvas->GetShowToolBar() ) cCanvas->ToggleToolBar(); + cCanvas->Modified(); cCanvas->Update(); + cCanvas->cd(); cCanvas->Divide(canvasDiv[0],canvasDiv[1]); + + gStyle->SetOptStat("neiou"); + + int padID = 0; + + for( int i = 0; i < NCLOVER; i++){ + padID++; + cCanvas->cd(padID); + cCanvas->cd(padID)->SetLogz(0); + hg[i]->Draw(""); + hg_g1[i]->SetLineColor(2);hg_g1[i]->Draw("same"); + hg_g2[i]->SetLineColor(4);hg_g2[i]->Draw("same"); + hg_g3[i]->SetLineColor(6);hg_g3[i]->Draw("same"); + } + + TCanvas *cPID = new TCanvas("cPID", "" ,7*200,4*200); + if( !cPID->GetShowToolBar() ) cPID->ToggleToolBar(); + cPID->Modified(); cPID->Update(); + cPID->cd(); cPID->Divide(7,4); + + for( int i = 0; i < NGAGG; i++){ + cPID->cd(i+1); + hPID[i]->Draw("colz"); + if( cutList1 != NULL && i < cutList1->GetEntries()) ((TCutG*)cutList1->At(i))->Draw("same"); + if( cutList2 != NULL && i < cutList2->GetEntries()) ((TCutG*)cutList2->At(i))->Draw("same"); + if( cutList3 != NULL && i < cutList3->GetEntries()) ((TCutG*)cutList3->At(i))->Draw("same"); + } + + + printf("=============== loaded AutoFit.C, try showFitMethos()\n"); + gROOT->ProcessLine(".L armory/AutoFit.C"); + + +} diff --git a/PIDAnalyzer.h b/PIDAnalyzer.h new file mode 100644 index 0000000..0013826 --- /dev/null +++ b/PIDAnalyzer.h @@ -0,0 +1,138 @@ +////////////////////////////////////////////////////////// +// This class has been automatically generated on +// Fri Mar 18 18:03:07 2022 by ROOT version 6.26/00 +// from TTree tree/tree +// found on file: haha.root +////////////////////////////////////////////////////////// + +#ifndef PIDAnalyzer_h +#define PIDAnalyzer_h + +#include +#include +#include +#include +#include + +#include "mapping.h" + +// Header file for the classes stored in the TTree if any. + +#define MAX_MULTI 40 + +class PIDAnalyzer : public TSelector { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + + // Declaration of leaf types + ULong64_t eventID; + Int_t runID; + Int_t multi; + Int_t multiGagg; + Short_t gammaID[MAX_MULTI]; //[multi] + Double_t gamma[MAX_MULTI]; //[multi] + ULong64_t gamma_t[MAX_MULTI]; //[multi] + Int_t gaggID[MAX_MULTI]; //[multiGagg] + Double_t gaggP[MAX_MULTI]; //[multiGagg] + Double_t gaggT[MAX_MULTI]; //[multiGagg] + ULong64_t gagg_t[MAX_MULTI]; //[multiGagg] + + // List of branches + TBranch *b_eventID; //! + TBranch *b_runID; //! + TBranch *b_multi; //! + TBranch *b_multiGagg; //! + TBranch *b_gammaID; //! + TBranch *b_gamma; //! + TBranch *b_gamma_t; //! + TBranch *b_gaggID; //! + TBranch *b_gaggP; //! + TBranch *b_gaggT; //! + TBranch *b_gagg_t; //! + + PIDAnalyzer(TTree * /*tree*/ =0) : fChain(0) { } + virtual ~PIDAnalyzer() { } + virtual Int_t Version() const { return 2; } + virtual void Begin(TTree *tree); + virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); + virtual Bool_t Notify(); + virtual Bool_t Process(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; } + virtual void SetOption(const char *option) { fOption = option; } + virtual void SetObject(TObject *obj) { fObject = obj; } + virtual void SetInputList(TList *input) { fInput = input; } + virtual TList *GetOutputList() const { return fOutput; } + virtual void SlaveTerminate(); + virtual void Terminate(); + + ClassDef(PIDAnalyzer,0); + + ULong64_t totnumEntry; + + ULong64_t ProcessedEntries; + Float_t Frac; ///Progress bar + TStopwatch StpWatch; +}; + +#endif + +#ifdef PIDAnalyzer_cxx +void PIDAnalyzer::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("eventID", &eventID, &b_eventID); + fChain->SetBranchAddress("runID", &runID, &b_runID); + fChain->SetBranchAddress("multi", &multi, &b_multi); + fChain->SetBranchAddress("multiGagg", &multiGagg, &b_multiGagg); + fChain->SetBranchAddress("gammaID", gammaID, &b_gammaID); + fChain->SetBranchAddress("gamma", gamma, &b_gamma); + fChain->SetBranchAddress("gamma_t", gamma_t, &b_gamma_t); + fChain->SetBranchAddress("gaggID", gaggID, &b_gaggID); + fChain->SetBranchAddress("gaggP", gaggP, &b_gaggP); + fChain->SetBranchAddress("gaggT", gaggT, &b_gaggT); + fChain->SetBranchAddress("gagg_t", gagg_t, &b_gagg_t); + + printf("======================== Start processing....\n"); + StpWatch.Start(); + + Frac = 0.1; + ProcessedEntries = 0; + +} + +Bool_t PIDAnalyzer::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + return kTRUE; +} + +void PIDAnalyzer::SlaveBegin(TTree * /*tree*/){ + + TString option = GetOption(); + +} + + +void PIDAnalyzer::SlaveTerminate(){ + +} + +#endif // #ifdef PIDAnalyzer_cxx diff --git a/PreAnalyzer.C b/PreAnalyzer.C new file mode 100644 index 0000000..af887b0 --- /dev/null +++ b/PreAnalyzer.C @@ -0,0 +1,195 @@ +#define PreAnalyzer_cxx + +#include "PreAnalyzer.h" +#include +#include + +#include + + +//############################################ BEGIN +void PreAnalyzer::Begin(TTree * tree){ + + TString option = GetOption(); + + totnumEntry = tree->GetEntries(); + + printf( "=========================================================================== \n"); + printf( "========================== PreAnalysis.C/h ============================= \n"); + printf( "====== total Entry : %lld \n", totnumEntry); + printf( "=========================================================================== \n"); + + +} + + + +//############################################ PROCESS +Bool_t PreAnalyzer::Process(Long64_t entry){ + + ProcessedEntries++; + + /*********** Progress Bar ******************************************/ + 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,totnumEntry/1000,StpWatch.RealTime(), StpWatch.RealTime()/Frac); + StpWatch.Start(kFALSE); + Frac+=0.1; + } + + b_event_ID->GetEntry(entry); + b_energy->GetEntry(entry); + b_time->GetEntry(entry); + b_multi->GetEntry(entry); + b_multiCry->GetEntry(entry); + b_detID->GetEntry(entry); + b_qdc->GetEntry(entry); + b_pileup->GetEntry(entry); + b_runID->GetEntry(entry); + + multi_N = 0; + multiGagg_N = 0; + + eventID = evID; + runID_N = runID; + + for( int i = 0; i < NCLOVER; i++) { + gammaID[i] = -1; + gamma_N[i] = TMath::QuietNaN(); + gamma_t[i] = 0; + } + + for( int i = 0; i < NGAGG; i++){ + gaggID[i] = -1; + gagg_peak[i] = TMath::QuietNaN(); + gagg_tail[i] = TMath::QuietNaN(); + gagg_t[i] = 0; + } + + double bg[NGAGG][2]={TMath::QuietNaN()}, peak[NGAGG][2]={TMath::QuietNaN()}, tail[NGAGG][2] = {TMath::QuietNaN()}; + ULong64_t gaggTime[NGAGG] = {0}; + ULong64_t gammaTime[NCLOVER] = {0}; + int count[NGAGG] = {0} ; + + ///################## Gamma data from Clover + for( int i = 0; i < NCRYSTAL; i++) eCal[i] = TMath::QuietNaN(); + for( int i = 0; i < NCLOVER; i++) gamma[i] = 0; + + for( int i = 0; i < multi ; i ++){ + if( pileup[i] == 1 ) continue; + int id = detID[i]; + + // GAGG_A + if( (200 <= id && id < 250) ) { + int id1 = id - 200; + + bg[id1][0] = (qdc[i][0] + qdc[i][1])/60.; + peak[id1][0] = qdc[i][3]/20. - bg[id1][0]; + tail[id1][0] = qdc[i][5]/55. - bg[id1][0]; + + if( gaggTime[id1] == 0 || e_t[i] < gaggTime[id1] ) gaggTime[id1] = e_t[i]; + + count[id1] ++; + } + + // GAGG_B + if( 250 <= id && id < 300 ) { + int id2 = id - 250; + + bg[id2][1] = (qdc[i][0] + qdc[i][1])/60.; + peak[id2][1] = qdc[i][3]/20. - bg[id2][1]; + tail[id2][1] = qdc[i][5]/55. - bg[id2][1]; + + if( gaggTime[id2] == 0 || e_t[i] < gaggTime[id2] ) gaggTime[id2] = e_t[i]; + + + count[id2]++; + } + + if( id >= 200 ) continue; + + //======== BGO veto + bool dropflag = false; + if( id < NCRYSTAL && multi > 1) { + for( int j = i + 1; j < multi; j++){ + if( 200 > detID[j] && detID[j] >= 100 && (detID[j]-100)*4 <= id && id < (detID[j]-100 +1)*4) { + dropflag = true; + break; + } + } + } + if( dropflag ) continue; + + if( 0<= id && id < NCRYSTAL ) { + if( eCorr.size() == 0 ){ + eCal[id] = e[i]; + }else{ + ///========= apply energy correction + int order = (int) eCorr[id].size(); + eCal[id] = 0; + for( int k = 0; k < order ; k++){ + eCal[id] += eCorr[id][k] * TMath::Power(e[i], k); + } + } + + ///========== add back + int cloverID = id /4; + if( eCal[id] > 10. ) { + gamma[cloverID] += eCal[id]; + if( gammaTime[cloverID] == 0 || e_t[i] < gammaTime[cloverID] ) gammaTime[cloverID] = e_t[i]; + + } + ///========= remove cross talk + + ///========= doppler correction + + } + + } + + + //################ Gamma-Paritcle + for( int i = 0 ; i < NCLOVER; i++){ + if( gamma[i] > 0 ) { + gammaID[multi_N] = i; + gamma_N[multi_N] = gamma[i]; + gamma_t[multi_N] = gammaTime[i]; + multi_N++; + } + } + + for ( int i = 0 ; i < NGAGG ; i++){ + if( count[i] == 2 ){ + gaggID[multiGagg_N] = i; + gagg_tail[multiGagg_N] = (tail[i][0]+tail[i][1])/2.; + gagg_peak[multiGagg_N] = (peak[i][0]+peak[i][1])/2.; + gagg_t[multiGagg_N] = gaggTime[i]; + multiGagg_N++; + } + } + + if( multi_N == 0 && multiGagg_N == 0 ) return kTRUE; + + //############################ save + saveFile->cd(); ///set focus on this file + newTree->Fill(); + + + return kTRUE; +} +//############################################ TERMINATE +void PreAnalyzer::Terminate(){ + + printf("============================== finishing.\n"); + + saveFile->cd(); //set focus on this file + newTree->Write(); + saveFile->Close(); + + printf("-------------- done, saved in %s.\n", saveFileName.Data()); + + gROOT->ProcessLine(".q"); + +} diff --git a/PreAnalyzer.h b/PreAnalyzer.h new file mode 100644 index 0000000..7330138 --- /dev/null +++ b/PreAnalyzer.h @@ -0,0 +1,200 @@ +////////////////////////////////////////////////////////// +// This class has been automatically generated on +// Mon Dec 13 18:37:46 2021 by ROOT version 6.24/06 +// from TTree tree/tree +// found on file: efEu152b-000.root +////////////////////////////////////////////////////////// + +#ifndef PreAnalyzer_h +#define PreAnalyzer_h + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "mapping.h" +#include "armory/AnalysisLibrary.h" + +// Header file for the classes stored in the TTree if any.' + +#define MAX_MULTI 200 + +class PreAnalyzer : public TSelector { +public : + TTree *fChain; //!pointer to the analyzed TTree or TChain + +// Fixed size dimensions of array or collections stored in the TTree if any. + + // Declaration of leaf types + ULong64_t evID; + Int_t runID; + 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; + Int_t multiGagg; + Bool_t pileup[MAX_MULTI]; + + // List of branches + TBranch *b_event_ID; //! + TBranch *b_runID; //! + TBranch *b_detID; //! + TBranch *b_energy; //! + TBranch *b_time; //! + TBranch *b_qdc; //! + TBranch *b_multi; //! + TBranch *b_multiCry; //! + TBranch *b_multiGagg; //! + TBranch *b_pileup; //! + + PreAnalyzer(TTree * /*tree*/ =0) : fChain(0) { totnumEntry = 0; + Frac = 0.1; + ProcessedEntries = 0; + } + virtual ~PreAnalyzer() { } + virtual Int_t Version() const { return 2; } + virtual void Begin(TTree *tree); + virtual void SlaveBegin(TTree *tree); + virtual void Init(TTree *tree); + virtual Bool_t Notify(); + virtual Bool_t Process(Long64_t entry); + virtual Int_t GetEntry(Long64_t entry, Int_t getall = 0) { return fChain ? fChain->GetTree()->GetEntry(entry, getall) : 0; } + virtual void SetOption(const char *option) { fOption = option; } + virtual void SetObject(TObject *obj) { fObject = obj; } + virtual void SetInputList(TList *input) { fInput = input; } + virtual TList *GetOutputList() const { return fOutput; } + virtual void SlaveTerminate(); + virtual void Terminate(); + + ClassDef(PreAnalyzer,0); + + ULong64_t totnumEntry; + + vector> eCorr; + + ULong64_t ProcessedEntries; + Float_t Frac; ///Progress bar + TStopwatch StpWatch; + + double eCal[NCRYSTAL]; + double gamma[NCLOVER]; // added-back energy + + ///========================= for new root file + TFile * saveFile; + TTree * newTree; + TString saveFileName; + + ///tree + ULong_t eventID; + Int_t runID_N; + int multi_N; + int multiGagg_N; + int gammaID[NCLOVER]; + double gamma_N[NCLOVER]; + ULong64_t gamma_t[NCLOVER]; + int gaggID[NGAGG]; + double gagg_peak[NGAGG]; + double gagg_tail[NGAGG]; + ULong64_t gagg_t[NGAGG]; + +}; + +#endif + +#ifdef PreAnalyzer_cxx +void PreAnalyzer::Init(TTree *tree) +{ + // The Init() function is called when the selector needs to initialize + // a new tree or chain. Typically here the branch addresses and branch + // pointers of the tree will be set. + // It is normally not necessary to make changes to the generated + // code, but the routine can be extended by the user if needed. + // Init() will be called many times when running on PROOF + // (once per file to be processed). + + // Set branch addresses and branch pointers + if (!tree) return; + fChain = tree; + fChain->SetMakeClass(1); + + fChain->SetBranchAddress("evID", &evID, &b_event_ID); + fChain->SetBranchAddress("runID", &runID, &b_runID); + 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); + fChain->SetBranchAddress("multiGagg", &multiGagg, &b_multiGagg); + fChain->SetBranchAddress("pileup", pileup, &b_pileup); + + TString option = GetOption(); + + printf("======================== Load parameters.\n"); + eCorr = LoadCorrectionParameters("correction_e.dat"); + + ///======================== open a new file + saveFileName = "haha.root"; + + saveFile = new TFile( saveFileName,"recreate"); + + TMacro e_corr("correction_e.dat"); + e_corr.Write("correction_e"); + + newTree = new TTree("tree", "tree"); + + eventID = -1; + runID = 0; + + multi_N = 0; + multiGagg_N = 0; + + newTree->Branch("eventID", &eventID, "eventID/l"); + newTree->Branch("runID", &runID_N, "runID/I"); + newTree->Branch("multi", &multi_N, "multi/I"); + newTree->Branch("multiGagg", &multiGagg_N, "multiGagg/I"); + newTree->Branch("gammaID", gammaID, "gammaID[multi]/S"); + newTree->Branch("gamma", gamma_N, "gamma[multi]/D"); + newTree->Branch("gamma_t", gamma_t, "gamma_t[multi]/l"); + newTree->Branch("gaggID", gaggID, "gaggID[multiGagg]/I"); + newTree->Branch("gaggP", gagg_peak, "gaggP[multiGagg]/D"); + newTree->Branch("gaggT", gagg_tail, "gaggT[multiGagg]/D"); + newTree->Branch("gagg_t", gagg_t, "gagg_t[multiGagg]/l"); + + printf("======================== Start processing....\n"); + StpWatch.Start(); + +} + +Bool_t PreAnalyzer::Notify() +{ + // The Notify() function is called when a new file is opened. This + // can be either for a new TTree in a TChain or when when a new TTree + // is started when using PROOF. It is normally not necessary to make changes + // to the generated code, but the routine can be extended by the + // user if needed. The return value is currently not used. + + return kTRUE; +} + + + +void PreAnalyzer::SlaveBegin(TTree * /*tree*/){ + + TString option = GetOption(); + +} + + +void PreAnalyzer::SlaveTerminate(){ + +} + +#endif // #ifdef PreAnalyzer_cxx