rtang/XIAEventBuilder
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
XIAEventBuilder/Analyzer.C

526 lines
17 KiB
C
Raw Blame History

#define Analyzer_cxx
#include "Analyzer.h"
#include <TH2.h>
#include <TStyle.h>
#include <TH1.h>
#include <TCutG.h>
#include <TCanvas.h>
#include <TMath.h>
#include <vector>
#include <stdio.h>
//############################################ User setting
int rawEnergyRange[2] = {0, 12000}; // in ch
int energyRange[3] = {1, 30, 800}; // keV {resol, min, max}
double BGO_threshold = 0; // in ch
int pidMaxRange[3] = {500, 400, 1600}; //nBin, tail, and peak
TString e_corr = "correction_e.dat";
TString cutFileName1 = "protonCut.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}
int runIDOffset = 1004000;
bool save_ev2 = false;
//############################################ end of user setting
//############################################ histogram declaration
TH2F * heVID;
TH1F * he[NCRYSTAL];
//TH2F * hgg[NCRYSTAL][NCRYSTAL];
TH2F * hcoin;
TH2F * hcrystalBGO;
TH2F * hgg;
TH1F * hTDiff;
TH2F * hPID[NGAGG];
TH2F * hPIDm[NGAGG];
TH2F * hPID_A[NGAGG];
TH2F * hPID_B[NGAGG];
TH2F * hGAGGVID;
///----- after gamma gate
TH2F * hPID_A_g[NGAGG];
TH2F * hpeakVrun;
TH2F * htailVrun;
///----- after calibration and BGO veto
TH2F * heCalVID;
TH1F * heCal[NCRYSTAL];
TH2F * hcoinBGO;
TH2F * hcrystalBGO_G;
TH1F * hg[NCLOVER];
///----- after particle gate
TH1F * heCal_g[NCRYSTAL];
TH2F * heCalVID_g;
TH1F * hg_g1[NCLOVER];
TH1F * hg_g2[NCLOVER];
TH1F * hg_g3[NCLOVER];
///============= cut
TCutG * cut1;
TCutG * cut2;
TCutG * cut3;
//############################################ BEGIN
void Analyzer::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");
heVID = new TH2F("heVID", "e vs ID; det ID; e [ch]", NCRYSTAL, 0, NCRYSTAL, rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
heCalVID = new TH2F("heCalVID", Form("eCal vs ID (BGO veto > %.1f); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
heCalVID_g = new TH2F("heCalVID_g", Form("eCal vs ID (BGO veto > %.1f & particle); det ID; Energy [keV]", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
hTDiff = new TH1F("hTDiff", "data time different within an event; tick [10 ns]", 110, 0, 110);
heVID->SetNdivisions(-410, "X");
heCalVID->SetNdivisions(-410, "X");
for( int i = 0; i < NCRYSTAL; i ++){
he[i] = new TH1F( Form("he%02d", i), Form("e -%02d", i), rawEnergyRange[1] - rawEnergyRange[0], rawEnergyRange[0], rawEnergyRange[1]);
heCal[i] = new TH1F(Form("heCal%02d", i), Form("eCal-%02d (BGO veto > %.1f); Energy [keV]; count / %d keV", i, BGO_threshold, energyRange[0]), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
heCal_g[i] = new TH1F(Form("heCal_g%02d", i), Form("eCal-%02d (BGO veto > %.1f & particle); Energy [keV]; count / %d keV", i, BGO_threshold, energyRange[0]), (energyRange[2] - energyRange[1])/energyRange[0], energyRange[1], energyRange[2]);
}
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]);
hPIDm[i] = new TH2F(Form("hPIDm%02d", i), Form("PIDm-%2d; full; peak - tail ", i), pidMaxRange[0], -20, pidMaxRange[1]+2000, pidMaxRange[0], -50, pidMaxRange[2]);
hPID_A[i] = new TH2F(Form("hPID_A%02d",i), Form("PID_A GAGG = %2d; tail; peak ", i) , pidMaxRange[0], -20, pidMaxRange[1], pidMaxRange[0], -50, pidMaxRange[2]);
hPID_B[i] = new TH2F(Form("hPID_B%02d",i), Form("PID_B GAGG = %2d; tail; peak ", i) , pidMaxRange[0], -20, pidMaxRange[1], pidMaxRange[0], -50, pidMaxRange[2]);
hPID_A_g[i] = new TH2F(Form("hPID_A_g%02d",i), Form("PID_A GAGG = %2d (gated); tail; peak ", i) , pidMaxRange[0], -20, pidMaxRange[1], pidMaxRange[0], -50, pidMaxRange[2]);
}
hpeakVrun = new TH2F("hpeakVrun", "GAGG-9 peak vs run", 200, 0, 200, pidMaxRange[0], -50, pidMaxRange[2]);
htailVrun = new TH2F("htailVrun", "GAGG-9 tail vs run", 200, 0, 200, pidMaxRange[0], -50, pidMaxRange[1]);
hGAGGVID = new TH2F("hGAGGVID", "GAGG V ID", 80, 0, 80, 400, -50, 2000);
/**
for( int i = 0; i < NCRYSTAL; i++){
for( int j = i+1; j < NCRYSTAL; j++){
//hgg[i][j] = new TH2F(Form("hgg%02d%02d", i, j), Form("e%02d vs e%02d; e%02d; e%02d", i, j, i, j),
// (rawEnergyRange[1] - rawEnergyRange[0])/2, rawEnergyRange[0], rawEnergyRange[1],
// (rawEnergyRange[1] - rawEnergyRange[0])/2, rawEnergyRange[0], rawEnergyRange[1]);
}
}*/
hcoin = new TH2F("hcoin", "detector coin.; det ID; det ID", NCRYSTAL, 0, NCRYSTAL, NCRYSTAL, 0 , NCRYSTAL);
hcoinBGO = new TH2F("hcoinBGO", Form("detector coin. (BGO veto > %.1f); det ID; det ID", BGO_threshold), NCRYSTAL, 0, NCRYSTAL, NCRYSTAL, 0 , NCRYSTAL);
hcrystalBGO = new TH2F("hcrystalBGO", Form("crystal vs BGO ; det ID; BGO ID"), NCRYSTAL, 0, NCRYSTAL, NBGO, 0 , NBGO);
hcrystalBGO_G = new TH2F("hcrystalBGO_G", Form("crystal vs BGO (BGO veto); det ID; BGO ID"), NCRYSTAL, 0, NCRYSTAL, NBGO, 0 , NBGO);
printf("======================== Load parameters.\n");
eCorr = LoadCorrectionParameters(e_corr);
printf("======================== Load cuts.\n");
if( cutFileName1 != ""){
TFile * cutFile1 = new TFile(cutFileName1);
cut1 = (TCutG *) cutFile1->Get("CUTG");
cut1->SetLineColor(2);
printf("%s| %s is loaded.\n", cutFileName1.Data(),cut1->GetName());
}
if( cutFileName2 != ""){
TFile * cutFile2 = new TFile(cutFileName2);
cut2 = (TCutG *) cutFile2->Get("CUTG");
cut2->SetLineColor(6);
printf("%s| %s is loaded.\n", cutFileName2.Data(),cut2->GetName());
}
if( cutFileName3 != ""){
TFile * cutFile3 = new TFile(cutFileName3);
cut3 = (TCutG *) cutFile3->Get("CUTG");
cut3->SetLineColor(5);
printf("%s| %s is loaded.\n", cutFileName3.Data(),cut3->GetName());
}
saveEV2 = save_ev2;
}
void Analyzer::PID_calculation(int ID){
}
//############################################ PROCESS
Bool_t Analyzer::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_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);
//################### Make PID from GAGG
double bg[NGAGG][2]={TMath::QuietNaN()}, peak[NGAGG][2]={TMath::QuietNaN()}, tail[NGAGG][2] = {TMath::QuietNaN()};
int count[NGAGG] = {0} ;
for( int i = 0; i < multi ; i ++){
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;
if( pileup[i] == 1 ) continue;
if( e[i] < 100 ) continue;
int id = detID[i];
if( id < 200 || id >= 300 ) continue;
id = id - 200;
hGAGGVID->Fill(id, e[i]);
// GAGG_A
if( (0 <= id && id < 50) ) {
bg[id][0] = (qdc[i][0] + qdc[i][1])/60.;
peak[id][0] = qdc[i][3]/20. - bg[id][0];
tail[id][0] = qdc[i][5]/55. - bg[id][0];
hPID_A[id]->Fill( tail[id][0], peak[id][0]);
count[id] ++;
}
// GAGG_B
if( 50 <= id ) {
id = id - 50;
bg[id][1] = (qdc[i][0] + qdc[i][1])/60.;
peak[id][1] = qdc[i][3]/20. - bg[id][1];
tail[id][1] = qdc[i][5]/55. - bg[id][1];
hPID_B[id]->Fill( tail[id][1], peak[id][1]);
count[id]++;
}
}
double tailAvg[NGAGG] = {TMath::QuietNaN()};
double peakAvg[NGAGG] = {TMath::QuietNaN()};
for ( int i = 0 ; i < NGAGG ; i++){
if( count[i] == 2 ){
tailAvg[i] = (tail[i][0]+tail[i][1])/2.;
peakAvg[i] = (peak[i][0]+peak[i][1])/2.;
hPID[i]->Fill( tailAvg[i], peakAvg[i]);
hPIDm[i]->Fill( peakAvg[i]-tailAvg[i], 100*(peakAvg[i] / tailAvg[i]) );
if( i == 9 ){
hpeakVrun->Fill( (runID - runIDOffset)/10 + (runID - runIDOffset)%1000, peakAvg[i] );
htailVrun->Fill( (runID - runIDOffset)/10 + (runID - runIDOffset)%1000, tailAvg[i] );
}
}
}
///################## 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;
if( 200 < detID[i] ) continue;
int id = detID[i];
//======== Fill raw data
if( 0 <= id && id < NCRYSTAL ){ /// gamma data
heVID->Fill( id, e[i]);
he[id]->Fill(e[i]);
for ( int j = i + 1; j < multi; j++){
if( 100 <= detID[j] && detID[j] < 200 ) hcrystalBGO->Fill(id, detID[j]-100); /// crystal - BGO coincident
if( detID[j] < 100 ) hcoin->Fill(id, detID[j]); /// crystal-crystal coincident
}
}
//======== 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;
if ( i > 0 ) hTDiff->Fill( e_t[i] - e_t[0]);
//======== BGO veto
bool dropflag = false;
if( id < NCRYSTAL && multi > 1) {
for( int j = i + 1; j < multi; j++){
if( 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( e_corr == "" ){
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);
}
}
if( id != 30 ) heCalVID->Fill( id, eCal[id]);
heCal[id]->Fill(eCal[id]);
for ( int j = i + 1; j < multi; j++){
if( 100 <= detID[j] && detID[j] < 200 ) hcrystalBGO_G->Fill(id, detID[j]-100); /// crystal - BGO coincident
}
///========== add back
int cloverID = id /4;
if( eCal[id] > energyRange[1]/4. ) gamma[cloverID] += eCal[id];
///========= remove cross talk
///========= doppler correction
}
}
//################ Gamma-Paritcle
bool fillFlag1 = false;
bool fillFlag2 = false;
bool fillFlag3 = false;
for( int gi = 0; gi < NGAGG ; gi ++){
if( fillFlag1 == false && cut1 != NULL && cut1->IsInside(tailAvg[gi], peakAvg[gi]) ) {
fillFlag1 = true;
}
if( fillFlag2 == false && cut2 != NULL && cut2->IsInside(tailAvg[gi], peakAvg[gi]) ) {
fillFlag2 = true;
}
if( fillFlag2 == false && cut3 != NULL && cut3->IsInside(tailAvg[gi], peakAvg[gi]) ) {
fillFlag3 = true;
}
}
for( int i = 0; i < NCLOVER; i++){
for( int j = 0; j < NCLOVER; j++){
if( gamma[i] > 0 && gamma[j] > 0 && i != j && fillFlag2) {
hgg->Fill( gamma[i], gamma[j]);
}
}
}
for( int i = 0 ; i < NCLOVER; i++){
if( gamma[i] > 0 ) {
hg[i]->Fill(gamma[i]);
if( fillFlag1 ) hg_g1[i]->Fill(gamma[i]);
if( fillFlag2 ) hg_g2[i]->Fill(gamma[i]);
if( fillFlag3 ) hg_g3[i]->Fill(gamma[i]);
}
if( abs(gamma[i] - 1052 ) < 8 ){
for( int i = 0; i < NGAGG ; i ++){
hPID_A_g[i]->Fill( tail[i][0], peak[i][0]);
}
}
}
if ( saveEV2) Save2ev2();
return kTRUE;
}
//############################################ TERMINATE
void Analyzer::Terminate(){
if(saveEV2) fclose(outEV2);
printf("============================== finishing.\n");
gROOT->cd();
int canvasXY[2] = {1600 , 800} ;// x, y
int canvasDiv[2] = {4,2};
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;
//========================= canvas 1
padID++;
cCanvas->cd(padID);
cCanvas->cd(padID)->SetLogz(0);
heCalVID->Draw("colz");
//========================= canvas 2
padID++;
cCanvas->cd(padID);
cCanvas->cd(padID)->SetLogz(1);
hGAGGVID->Draw("colz");
//========================= canvas 3
padID++;
cCanvas->cd(padID);
cCanvas->cd(padID)->SetLogz(0);
hTDiff->Draw();
//========================= canvas 4
padID++;
cCanvas->cd(padID);
cCanvas->cd(padID)->SetLogz(0);
hg[6]->Draw();
//========================= canvas 5
padID++;
cCanvas->cd(padID);
cCanvas->cd(padID)->SetLogz(1);
//heCalVID->Draw("colz");
hPID[9]->Draw("colz");
if( cut1 != NULL ) cut1->Draw("same");
if( cut2 != NULL ) cut2->Draw("same");
//========================= canvas 6
padID++;
cCanvas->cd(padID);
cCanvas->cd(padID)->SetLogz(1);
//heCalVID->Draw("colz");
//hPIDm[9]->Draw("colz");
hpeakVrun->Draw("colz");
//hPID_A[9]->Draw("colz");
//========================= canvas 7
padID++;
cCanvas->cd(padID);
cCanvas->cd(padID)->SetLogz(1);
//htailVrun->Draw("colz");
hgg->Draw("colz");
//hPID_B[9]->Draw("colz");
//========================= canvas 8
padID++;
cCanvas->cd(padID);
cCanvas->cd(padID)->SetLogz(1);
//hPID_A_g[9]->Draw("colz");
hg_g2[6]->SetLineColor(2);
hg_g2[6]->Draw();
hg_g3[6]->SetLineColor(4);
hg_g3[6]->Draw("same");
hg_g1[6]->Draw("same");
/*
//========================= canvas 1
cCanvas->cd(4);
cCanvas->cd(4)->SetLogy(1);
//gROOT->ProcessLine(".x script.C");
//hcrystalBGO_G->Draw("colz");
hg[0]->SetLineColor(2);
hg[0]->Draw();
hg_g[0]->Draw("same");
//hcoinBGO->Draw("colz");
TCanvas *cAux = new TCanvas("cAux", "" ,1000, 0, 800,600);
cAux->cd();
TH1F * h0 = (TH1F *) heCal[0]->Clone("h0");
h0->Add(heCal[1]);
h0->Add(heCal[2]);
h0->Add(heCal[3]);
TH1F * h0_g = (TH1F *) heCal_g[0]->Clone("h0_g");
h0_g->Add(heCal_g[1]);
h0_g->Add(heCal_g[2]);
h0_g->Add(heCal_g[3]);
h0->SetLineColor(kGreen+3);
h0->Draw();
h0_g->Draw("same");
hg[0]->Draw("same");
*/
printf("=============== loaded AutoFit.C, try showFitMethos()\n");
gROOT->ProcessLine(".L armory/AutoFit.C");
printf("=============== Analyzer Utility\n");
gROOT->ProcessLine(".L armory/Analyzer_Utili.c");
gROOT->ProcessLine("listDraws()");
}
Reference in New Issue View Git Blame Copy Permalink