added GAGGPIDCutCreator, PreAnalyzer to reduce the root size, PIDAnalyzer
This commit is contained in:
nrb@Debain10
parent
3066de5c09
commit
9700accb0f
97
GAGGPIDCutCreator.C
Normal file
97
GAGGPIDCutCreator.C
Normal file
|
|
@ -0,0 +1,97 @@
|
|||
#include <TH2F.h>
|
||||
#include <TFile.h>
|
||||
#include <TChain.h>
|
||||
#include <TROOT.h>
|
||||
#include <TStyle.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TCutG.h>
|
||||
#include <TString.h>
|
||||
#include <TObjArray.h>
|
||||
|
||||
#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");
|
||||
|
||||
}
|
||||
263
PIDAnalyzer.C
Normal file
263
PIDAnalyzer.C
Normal file
|
|
@ -0,0 +1,263 @@
|
|||
#define PIDAnalyzer_cxx
|
||||
|
||||
|
||||
#include "PIDAnalyzer.h"
|
||||
#include <TH2.h>
|
||||
#include <TStyle.h>
|
||||
#include <TH1.h>
|
||||
#include <TCutG.h>
|
||||
#include <TCanvas.h>
|
||||
#include <TMath.h>
|
||||
#include <TObjArray.h>
|
||||
#include <vector>
|
||||
#include <stdio.h>
|
||||
|
||||
//############################################ 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");
|
||||
|
||||
|
||||
}
|
||||
138
PIDAnalyzer.h
Normal file
138
PIDAnalyzer.h
Normal file
|
|
@ -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 <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
#include <TStopwatch.h>
|
||||
|
||||
#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
|
||||
195
PreAnalyzer.C
Normal file
195
PreAnalyzer.C
Normal file
|
|
@ -0,0 +1,195 @@
|
|||
#define PreAnalyzer_cxx
|
||||
|
||||
#include "PreAnalyzer.h"
|
||||
#include <TStyle.h>
|
||||
#include <TMath.h>
|
||||
|
||||
#include <stdio.h>
|
||||
|
||||
|
||||
//############################################ 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");
|
||||
|
||||
}
|
||||
200
PreAnalyzer.h
Normal file
200
PreAnalyzer.h
Normal file
|
|
@ -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 <fstream>
|
||||
#include <TROOT.h>
|
||||
#include <TChain.h>
|
||||
#include <TFile.h>
|
||||
#include <TSelector.h>
|
||||
#include <TStopwatch.h>
|
||||
#include <TMacro.h>
|
||||
#include <vector>
|
||||
|
||||
#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<vector<double>> 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
|
||||
Loading…
Reference in New Issue
Block a user