268 lines
6.6 KiB
C
268 lines
6.6 KiB
C
#include <TFile.h>
|
|
#include <TTree.h>
|
|
#include <TCanvas.h>
|
|
#include <TROOT.h>
|
|
#include <TSystem.h>
|
|
#include <TStyle.h>
|
|
#include <TProfile.h>
|
|
#include <TH2F.h>
|
|
#include <TH1F.h>
|
|
#include <TF1.h>
|
|
#include <TMath.h>
|
|
#include <TGraph.h>
|
|
#include <TLine.h>
|
|
#include <TSpectrum.h>
|
|
#include "armory/AnalysisLibrary.h"
|
|
|
|
vector<double> Cali_gamma_histogram(TH1F * hist, int opt = 0, float threshold = 0.1, int peakDensity = 10){
|
|
/**///======================================================== load tree
|
|
|
|
printf("============================================================= \n");
|
|
printf("====================== Cali_gamma.C ========================= \n");
|
|
printf("============================================================= \n");
|
|
|
|
/**///======================================================== Browser or Canvas
|
|
|
|
vector<double> sol;
|
|
|
|
if( hist->GetEntries() < 100 ){
|
|
sol.push_back(0);
|
|
sol.push_back(1);
|
|
return sol;
|
|
}
|
|
|
|
Int_t Div[2] = {2,2}; //x,y
|
|
Int_t size[2] = {600,300}; //x,y
|
|
|
|
TCanvas * cAlpha = (TCanvas *) gROOT->FindObjectAny("cAlpha");
|
|
if( cAlpha == NULL ){
|
|
cAlpha = new TCanvas("cAlpha", "cAlpha", 0, 0, size[0]*Div[0], size[1]*Div[1]);
|
|
}
|
|
cAlpha->Clear();
|
|
cAlpha->Divide(Div[0],Div[1]);
|
|
|
|
for( int i = 1; i <= Div[0]*Div[1] ; i++){
|
|
cAlpha->cd(i)->SetGrid();
|
|
}
|
|
|
|
gStyle->SetOptStat(0);
|
|
gStyle->SetStatY(1.0);
|
|
gStyle->SetStatX(0.99);
|
|
gStyle->SetStatW(0.2);
|
|
gStyle->SetStatH(0.1);
|
|
|
|
if(cAlpha->GetShowEditor() )cAlpha->ToggleEditor();
|
|
if(cAlpha->GetShowToolBar() )cAlpha->ToggleToolBar();
|
|
|
|
/**///========================================================= Analysis
|
|
|
|
cAlpha->cd(1);
|
|
hist->Draw();
|
|
cAlpha->Update();
|
|
gSystem->ProcessEvents();
|
|
|
|
int dummy = 0;
|
|
int temp = 0;
|
|
|
|
int nPeaks = 10;
|
|
vector<double> energy;
|
|
vector<double> refEnergy;
|
|
vector<double> fitEnergy;
|
|
|
|
|
|
printf("---- finding peak using TSepctrum Class...\n");
|
|
|
|
TSpectrum * spec = new TSpectrum(20);
|
|
nPeaks = spec->Search(hist, peakDensity, "", threshold);
|
|
printf("---- found %2d peaks | ", nPeaks);
|
|
|
|
double * xpos = spec->GetPositionX();
|
|
double * ypos = spec->GetPositionY();
|
|
|
|
cAlpha->Update();
|
|
gSystem->ProcessEvents();
|
|
|
|
vector<double> height;
|
|
|
|
int * inX = new int[nPeaks];
|
|
TMath::Sort(nPeaks, xpos, inX, 0 );
|
|
for( int j = 0; j < nPeaks; j++){
|
|
energy.push_back(xpos[inX[j]]);
|
|
height.push_back(ypos[inX[j]]);
|
|
}
|
|
|
|
for( int j = 0; j < nPeaks; j++) printf("%7.2f, ", energy[j]);
|
|
printf("\n");
|
|
|
|
|
|
//------------ 3, correction
|
|
int option = 0;
|
|
if ( opt == 0 ) {
|
|
printf("========== which detector to be the reference?\n");
|
|
printf("-1 = manual reference\n");
|
|
printf("-2 = use 228Th, first 5 strongest peaks \n");
|
|
printf("-3 = use 207Bi, 3 peaks \n");
|
|
printf("-4 = use 152Eu, 11 peaks \n");
|
|
printf("-9 = stop \n");
|
|
printf("your choice = ");
|
|
temp = scanf("%d", &option);
|
|
}else{
|
|
option = opt;
|
|
}
|
|
|
|
if( option == -9 ) {
|
|
printf("------ stopped by user.\n");
|
|
return sol;
|
|
}
|
|
|
|
|
|
|
|
//======== fill reference energy
|
|
if(option == -1){
|
|
for( int i = 0; i < nPeaks; i++){
|
|
float eng;
|
|
printf("%2d-th peak energy %f ch (-1 to skip, -2 to skip all):", i, energy[i]);
|
|
temp = scanf("%f", &eng);
|
|
if( eng >= 0 ) {
|
|
refEnergy.push_back(eng);
|
|
fitEnergy.push_back(energy[i]);
|
|
printf(" input: %f \n", eng);
|
|
}else if ( eng == -1 ){
|
|
printf(" input: skipped \n");
|
|
}else if ( eng == -2 ){
|
|
break;
|
|
}
|
|
};
|
|
}
|
|
|
|
if( option == -2 ){
|
|
fitEnergy = energy;
|
|
|
|
refEnergy.clear();
|
|
refEnergy.push_back(5.423);
|
|
refEnergy.push_back(5.685);
|
|
refEnergy.push_back(6.288);
|
|
refEnergy.push_back(6.778);
|
|
refEnergy.push_back(8.785);
|
|
}
|
|
|
|
if( option == -3 ){
|
|
|
|
fitEnergy = energy;
|
|
|
|
refEnergy.clear();
|
|
refEnergy.push_back(0.569702);
|
|
refEnergy.push_back(1.063662);
|
|
refEnergy.push_back(1.770237);
|
|
}
|
|
|
|
if( option == -4 ){
|
|
fitEnergy = energy;
|
|
|
|
refEnergy.clear();
|
|
refEnergy.push_back( 121.783);
|
|
refEnergy.push_back( 244.699);
|
|
refEnergy.push_back( 344.281);
|
|
refEnergy.push_back( 411.115);
|
|
refEnergy.push_back( 443.965);
|
|
refEnergy.push_back( 778.903);
|
|
refEnergy.push_back( 867.390);
|
|
refEnergy.push_back( 964.055);
|
|
refEnergy.push_back(1085.842);
|
|
refEnergy.push_back(1112.087);
|
|
refEnergy.push_back(1408.022);
|
|
}
|
|
|
|
//==================== adjusting energy
|
|
int n = refEnergy.size();
|
|
for( int k = 0; k < n; k++) printf("%2d-th peak : %f \n", k, refEnergy[k]);
|
|
|
|
nPeaks = energy.size();
|
|
|
|
vector<double> k1 = fitEnergy;
|
|
vector<double> k2 = refEnergy;
|
|
//===== when nPeaks != refEnergy.size(), need to matching the two vector size by checking the r-squared.
|
|
if( option != -1 ){
|
|
vector<vector<double>> haha = FindMatchingPair(fitEnergy, refEnergy);
|
|
k1 = haha[0];
|
|
k2 = haha[1];
|
|
}
|
|
|
|
TGraph * graph = new TGraph(min(n, nPeaks), &k1[0], &k2[0] );
|
|
cAlpha->cd(2);
|
|
graph->Draw("A*");
|
|
gSystem->ProcessEvents();
|
|
|
|
TF1 * fit = new TF1("fit", "pol1" );
|
|
graph->Fit("fit", "q");
|
|
|
|
double a0 = fit->GetParameter(0);
|
|
double a1 = fit->GetParameter(1);
|
|
double a2 = 0.0 ; //fit->GetParameter(2);
|
|
|
|
printf("----- a0: %9.6f, a1: %9.6f (%14.8f) \n", a0, a1, 1./a1);
|
|
printf(" %13.10f\t %13.10f \n", a0, a1);
|
|
|
|
sol.clear();
|
|
sol.push_back(a0);
|
|
sol.push_back(a1);
|
|
|
|
//====== Plot residue
|
|
vector<double> resi;
|
|
for( int i = 0 ; i < min(n, nPeaks); i++){
|
|
resi.push_back(k2[i] - fit->Eval(k1[i]));
|
|
}
|
|
TGraph * graphResi = new TGraph(min(n, nPeaks), &k2[0], &resi[0] );
|
|
|
|
cAlpha->cd(4);
|
|
graphResi->Draw("A*");
|
|
graphResi->GetYaxis()->SetTitle("residue [keV]");
|
|
graphResi->GetXaxis()->SetTitle("cali energy [keV]");
|
|
gSystem->ProcessEvents();
|
|
|
|
//====== Plot adjusted spectrum
|
|
|
|
cAlpha->cd(3);
|
|
|
|
int bin = hist->GetNbinsX();
|
|
double xMin = hist->GetXaxis()->GetXmin();
|
|
double xMax = hist->GetXaxis()->GetXmax();
|
|
|
|
double xMinC, xMaxC;
|
|
if( a1 > 0 ) {
|
|
xMinC = xMin*xMin*a2 + xMin*a1 + a0;
|
|
xMaxC = xMax*xMax*a2 + xMax*a1 + a0;
|
|
}else{
|
|
xMaxC = xMin*xMin*a2 + xMin*a1 + a0;
|
|
xMinC = xMax*xMax*a2 + xMax*a1 + a0;
|
|
}
|
|
|
|
|
|
TH1F * calH = new TH1F("calH", "calibrated energy", bin, xMinC, xMaxC);
|
|
|
|
/*
|
|
FILE * paraOut;
|
|
TString filename;
|
|
filename.Form("hist%s.dat", hist->GetName());
|
|
paraOut = fopen (filename.Data(), "w+");
|
|
*/
|
|
for( int i = 1; i <= bin; i ++){
|
|
int y = hist->GetBinContent(i);
|
|
int x = hist->GetBinCenter(i);
|
|
|
|
calH->Fill(x*x*a2 + x*a1+a0, y);
|
|
|
|
//fprintf(paraOut, "%9.6f, %9d\n", x*x*a2 + x*a1+a0, y);
|
|
|
|
}
|
|
|
|
// fflush(paraOut);
|
|
// fclose(paraOut);
|
|
|
|
calH->Draw();
|
|
gSystem->ProcessEvents();
|
|
|
|
return sol;
|
|
|
|
}
|