249 lines
7.6 KiB
C
249 lines
7.6 KiB
C
#define Analyzer_cxx
|
|
|
|
#include "Analyzer.h"
|
|
#include <TH2.h>
|
|
#include <TStyle.h>
|
|
#include <TH1.h>
|
|
#include <TCanvas.h>
|
|
#include <TMath.h>
|
|
#include <vector>
|
|
#include <stdio.h>
|
|
|
|
//############################################ User setting
|
|
|
|
int rawEnergyRange[2] = {0, 6000}; // in ch
|
|
int energyRange[3] = {1, 0, 6000}; // keV {resol, min, max}
|
|
|
|
double BGO_threshold = 0; // in ch
|
|
|
|
TString e_corr = "";//"correction_e.dat";
|
|
|
|
//############################################ end of user setting
|
|
|
|
|
|
//############################################ histogram declaration
|
|
|
|
TH2F * hevID;
|
|
TH1F * he[NCRYSTAL];
|
|
|
|
TH2F * hgg[NCRYSTAL][NCRYSTAL];
|
|
|
|
TH2F * hcoin;
|
|
|
|
///----- after calibration and BGO veto
|
|
TH2F * heCalvID;
|
|
TH1F * heCal[NCRYSTAL];
|
|
TH2F * hcoinBGO;
|
|
|
|
TH2F * hcrystalBGO;
|
|
|
|
//############################################ 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]);
|
|
|
|
hevID->SetNdivisions(-409, "X");
|
|
heCalvID->SetNdivisions(-409, "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]);
|
|
}
|
|
|
|
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);
|
|
|
|
printf("======================== End of histograms declaration\n");
|
|
|
|
printf("======================== Load parameters.\n");
|
|
|
|
eCorr = LoadCorrectionParameters(e_corr);
|
|
|
|
}
|
|
//############################################ 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_pileup->GetEntry(entry);
|
|
b_hit->GetEntry(entry);
|
|
b_bgo->GetEntry(entry);
|
|
b_bgoTime->GetEntry(entry);
|
|
b_other->GetEntry(entry);
|
|
b_multi->GetEntry(entry);
|
|
|
|
if( multi == 0 ) return kTRUE;
|
|
|
|
int numGatedData=0;
|
|
vector<int> gatedID;
|
|
gatedID.clear();
|
|
double eCal[NCRYSTAL] ={0.0};
|
|
|
|
///=========== 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;
|
|
|
|
//======== Fill raw data
|
|
hevID->Fill( detID, e[detID]);
|
|
he[detID]->Fill(e[detID]);
|
|
|
|
for( int kk = 0 ; kk < NBGO; kk++){
|
|
if( bgo[kk] > 0 ) hcrystalBGO->Fill(detID, kk);
|
|
}
|
|
|
|
|
|
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;
|
|
}
|
|
}
|
|
|
|
//----- for ev2 file
|
|
if( saveEV2 ){
|
|
numGatedData ++;
|
|
gatedID.push_back(detID);
|
|
}
|
|
|
|
//========= 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);
|
|
}
|
|
|
|
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);
|
|
}
|
|
}
|
|
|
|
|
|
if ( saveEV2){
|
|
|
|
short * out0 = new short[1];
|
|
short * outa = new short[1];
|
|
short * outb = new short[1];
|
|
|
|
out0[0] = numGatedData;
|
|
fwrite(out0, 1, 1, outEV2);
|
|
|
|
for( int i = 0; i < (int) gatedID.size(); i++){
|
|
int id = gatedID[i];
|
|
outa[0] = id;
|
|
fwrite(outa, 1, 1, outEV2);
|
|
outb[0] = eCal[id];
|
|
fwrite(outb, 2, 1, outEV2);
|
|
}
|
|
|
|
fwrite(out0, 1, 1, outEV2);
|
|
|
|
/**
|
|
int len = (int) gatedID.size();
|
|
char out[2*len+2];
|
|
out[0] = numGatedData;
|
|
for( int i = 0; i < (int) gatedID.size(); i++){
|
|
int id = gatedID[i];
|
|
out[2*i+1] = id;
|
|
out[2*i+2] = eCal[id];
|
|
}
|
|
out[2*len+1]=numGatedData;
|
|
fwrite(out,3*len+2,sizeof(out),outEV2);
|
|
*/
|
|
}
|
|
|
|
|
|
return kTRUE;
|
|
}
|
|
//############################################ TERMINATE
|
|
void Analyzer::Terminate(){
|
|
|
|
if(saveEV2) fclose(outEV2);
|
|
|
|
printf("============================== finishing.\n");
|
|
gROOT->cd();
|
|
|
|
int canvasXY[2] = {1200 , 1200} ;// x, y
|
|
int canvasDiv[2] = {2,2};
|
|
TCanvas *cCanvas = new TCanvas("cCanvas", "" ,canvasXY[0],canvasXY[1]);
|
|
cCanvas->Modified(); cCanvas->Update();
|
|
cCanvas->cd(); cCanvas->Divide(canvasDiv[0],canvasDiv[1]);
|
|
|
|
gStyle->SetOptStat("neiou");
|
|
|
|
cCanvas->cd(1);
|
|
cCanvas->cd(1)->SetLogz(1);
|
|
hevID->Draw("colz");
|
|
|
|
cCanvas->cd(2);
|
|
cCanvas->cd(2)->SetLogz(1);
|
|
heCalvID->Draw("colz");
|
|
|
|
cCanvas->cd(3);
|
|
cCanvas->cd(3)->SetLogz(1);
|
|
hcrystalBGO->Draw("colz");
|
|
|
|
cCanvas->cd(4);
|
|
//cCanvas->cd(4)->SetLogz(1);
|
|
he[0]->SetLineColor(2);
|
|
he[0]->Draw();
|
|
heCal[0]->Draw("same");
|
|
//hcoinBGO->Draw("colz");
|
|
|
|
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()");
|
|
|
|
|
|
|
|
|
|
}
|