Skip to content
Snippets Groups Projects
Commit c759a156 authored by gypaos's avatar gypaos
Browse files

Add L/M branching ratio in fitting procedure

parent 547207a9
No related branches found
No related tags found
No related merge requests found
Pipeline #104151 passed
Stage: build-NPLib
Stage: build-NPSimulation
Stage: test
Hide whitespace changes
Inline Side-by-side
Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C
+ 85
− 31
View file @ c759a156
...@@ -31,6 +31,7 @@ ...@@ -31,6 +31,7 @@
#include "TLine.h" #include "TLine.h"
#include "TError.h" #include "TError.h"
#include "TROOT.h" #include "TROOT.h"
#include "TMath.h"
// C++ headers // C++ headers
#include <iostream> #include <iostream>
...@@ -54,6 +55,7 @@ static std::vector<Double_t> peakListForFit; ...@@ -54,6 +55,7 @@ static std::vector<Double_t> peakListForFit;
// functions // functions
void AddPeak(Double_t, Double_t, Double_t); void AddPeak(Double_t, Double_t, Double_t);
Double_t fpeaks(Double_t*, Double_t*); Double_t fpeaks(Double_t*, Double_t*);
Double_t fpeaks2(Double_t*, Double_t*);
Double_t gaussianPeak(Double_t*, Double_t*); Double_t gaussianPeak(Double_t*, Double_t*);
...@@ -89,7 +91,7 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -89,7 +91,7 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
// open output pdf file // open output pdf file
string label = sname.substr(4, 6); string label = sname.substr(4, 6);
label += sname.substr(16, 20); label += sname.substr(16, 21);
label += (isPside) ? "_pside" : "_nside"; label += (isPside) ? "_pside" : "_nside";
TString pdfname = Form("Analyse207Bi_%s.pdf", label.c_str()); TString pdfname = Form("Analyse207Bi_%s.pdf", label.c_str());
can->Print(Form("%s[", pdfname.Data())); can->Print(Form("%s[", pdfname.Data()));
...@@ -101,7 +103,8 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -101,7 +103,8 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
// define some functions // define some functions
auto fcalib = new TF1("fcalib", "pol1", 1024); auto fcalibFull = new TF1("fcalibFull", "pol1", 1024);
auto fcalibLocal = new TF1("fcalibLocal", "pol1", 1024);
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
// define vectors for calibration // define vectors for calibration
...@@ -142,7 +145,7 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -142,7 +145,7 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
for (Int_t n = 0; n < NCHANNELS; ++n) { // loop on channels for (Int_t n = 0; n < NCHANNELS; ++n) { // loop on channels
// treat all channels or specified one // treat all channels or specified one
if (channel < 0 || channel == n) { if (channel < 0 || channel == n) {
cout << "Analysing channel " << n << "\n"; cout << "\n--------------- Analysing channel " << n << " ---------------\n";
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
// get histogram // get histogram
...@@ -168,8 +171,8 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -168,8 +171,8 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
auto s = new TSpectrum(); auto s = new TSpectrum();
Int_t npeaks = s->Search(h, 6, "", 0.6); Int_t npeaks = s->Search(h, 6, "", 0.6);
h->DrawCopy(); h->DrawCopy();
if (npeaks == 1) cout << "pedestal found!\n"; if (npeaks == 1) cout << ". Pedestal found -> OK\n";
else cout << "PROBLEM with pedestal!\n"; else cout << ". Pedestal PROBLEM\n";
Double_t *xpeaks = s->GetPositionX(); Double_t *xpeaks = s->GetPositionX();
// fit // fit
h->Fit("gaus", "RQ", "", xpeaks[0]-20, xpeaks[0]+20); h->Fit("gaus", "RQ", "", xpeaks[0]-20, xpeaks[0]+20);
...@@ -200,7 +203,7 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -200,7 +203,7 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
// search peaks after background subtraction // search peaks after background subtraction
npeaks = s->Search(h, 6, "", 0.3); npeaks = s->Search(h, 6, "", 0.3);
h->DrawCopy(); h->DrawCopy();
cout << "\t-> found " << npeaks << " peaks"; cout << ". " << npeaks << " main transitions found ";
if (npeaks == 2) cout << "\t-> OK!\n"; if (npeaks == 2) cout << "\t-> OK!\n";
else cout << "\t-> PROBLEM!\n"; else cout << "\t-> PROBLEM!\n";
// get peak position // get peak position
...@@ -214,13 +217,13 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -214,13 +217,13 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
// try to guess which peak is detected and add other one // try to guess which peak is detected and add other one
// uses calibration parameters from previous strip calibration // uses calibration parameters from previous strip calibration
if (npeaks == 1) { if (npeaks == 1) {
Double_t energy = fcalib->GetParameter(0) + fcalib->GetParameter(1)*xpeaks[0]; Double_t energy = fcalibFull->GetParameter(0) + fcalibFull->GetParameter(1)*xpeaks[0];
Double_t energyGuess = 0; Double_t energyGuess = 0;
if (energy > mainLineEnergy[0]+200) { if (energy > mainLineEnergy[0]+200) {
energyGuess = (mainLineEnergy[0]-fcalib->GetParameter(0))/fcalib->GetParameter(1);; energyGuess = (mainLineEnergy[0]-fcalibFull->GetParameter(0))/fcalibFull->GetParameter(1);;
} }
else { else {
energyGuess = (mainLineEnergy[1]-fcalib->GetParameter(0))/fcalib->GetParameter(1);; energyGuess = (mainLineEnergy[1]-fcalibFull->GetParameter(0))/fcalibFull->GetParameter(1);;
} }
AddPeak(energyGuess, BACKGROUND_MIN, BACKGROUND_MAX); AddPeak(energyGuess, BACKGROUND_MIN, BACKGROUND_MAX);
} }
...@@ -230,9 +233,9 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -230,9 +233,9 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
// rough energy calibration // rough energy calibration
auto grcalib = new TGraph(peakList.size(), &peakList[0], &mainLineEnergy[0]); auto grcalib = new TGraph(peakList.size(), &peakList[0], &mainLineEnergy[0]);
grcalib->Fit("fcalib", "Q0"); grcalib->Fit("fcalibFull", "Q0");
grOffset->SetPoint(n, n, fcalib->GetParameter(0)); grOffset->SetPoint(n, n, fcalibFull->GetParameter(0));
grGain ->SetPoint(n, n, fcalib->GetParameter(1)*10); grGain ->SetPoint(n, n, fcalibFull->GetParameter(1)*10);
// loop on "480"- and "975"-keV features // loop on "480"- and "975"-keV features
// p = 0 -> 480 keV // p = 0 -> 480 keV
...@@ -242,15 +245,22 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -242,15 +245,22 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
can->cd(3+p); can->cd(3+p);
pad = (TPad*) can->FindObject(Form("can_%d", 3+p)); pad = (TPad*) can->FindObject(Form("can_%d", 3+p));
// determine range for fitting based on rough calibration // determine range for fitting based on rough calibration
Double_t cmin = (mainLineEnergy[p]-40 -fcalib->GetParameter(0)) / fcalib->GetParameter(1); Double_t cmin = (mainLineEnergy[p]-40 -fcalibFull->GetParameter(0)) / fcalibFull->GetParameter(1);
Double_t cmax = (mainLineEnergy[p]+120-fcalib->GetParameter(0)) / fcalib->GetParameter(1); Double_t cmax = (mainLineEnergy[p]+130-fcalibFull->GetParameter(0)) / fcalibFull->GetParameter(1);
std::cout << cmin << "\t" << cmax << "\n"; if (cmax > 1023) cmax = 1023; // prevents detection of spurious peak at channel 1024
// std::cout << cmin << "\t" << cmax << "\n";
h->GetXaxis()->SetRangeUser(cmin, cmax); h->GetXaxis()->SetRangeUser(cmin, cmax);
h->DrawCopy(); h->DrawCopy();
// search satellite peaks // search satellite peaks
npeaks = s->Search(h, 4, "", 0.1); npeaks = s->Search(h, 4, "", 0.1);
cout << "\t-> found " << npeaks << " peaks"; if (p == 0) {
cout << ". Low energy transitions\n";
}
else {
cout << ". High energy transitions\n";
}
cout << "\t. found " << npeaks << " peaks";
if (npeaks == 2) cout << "\t-> OK!\n"; if (npeaks == 2) cout << "\t-> OK!\n";
else cout << "\t-> PROBLEM!\n"; else cout << "\t-> PROBLEM!\n";
// get peak position // get peak position
...@@ -266,14 +276,18 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -266,14 +276,18 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
// where a "replica" peak close to the main peak is observed, and // where a "replica" peak close to the main peak is observed, and
// should be removed // should be removed
if (npeaks == 3) peakListForFit.erase(peakListForFit.begin()+1); if (npeaks == 3) peakListForFit.erase(peakListForFit.begin()+1);
// add CE-M component // add CE-M component based on local linear fit
Double_t channel = (MLineEnergy[p] -fcalib->GetParameter(0)) / fcalib->GetParameter(1); vector<double> energyLocal = {mainLineEnergy[p], LLineEnergy[p]};
auto grcalibLocal = new TGraph(peakListForFit.size(), &peakListForFit[0], &energyLocal[0]);
grcalibLocal->Fit("fcalibLocal", "Q0");
Double_t channel = (MLineEnergy[p] -fcalibLocal->GetParameter(0)) / fcalibLocal->GetParameter(1);
// Double_t channel = (MLineEnergy[p] -fcalibFull->GetParameter(0)) / fcalibFull->GetParameter(1);
AddPeak(channel, cmin, cmax); AddPeak(channel, cmin, cmax);
// define parameter list for fit // define parameter list for fit
std::vector<Double_t> paramList; std::vector<Double_t> paramList;
// width (sigma); common for all states // width (sigma); common for all states
paramList.push_back(4); paramList.push_back(6);
// position and amplitude // position and amplitude
for (UInt_t i = 0; i < peakListForFit.size(); ++i) { for (UInt_t i = 0; i < peakListForFit.size(); ++i) {
paramList.push_back(peakListForFit[i]); paramList.push_back(peakListForFit[i]);
...@@ -284,20 +298,26 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -284,20 +298,26 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
} }
// define fit function // define fit function
auto fit = new TF1("fit", fpeaks, cmin, cmax, paramList.size()); auto fit = new TF1("fit", fpeaks2, cmin, cmax, paramList.size());
fit->SetParameters(&paramList[0]); fit->SetParameters(&paramList[0]);
// positive amplitudes and position in range // positive amplitudes and position in range
int dChannel = 5;
for (UInt_t i = 0; i < peakListForFit.size(); ++i) { for (UInt_t i = 0; i < peakListForFit.size(); ++i) {
// +/- 10 channels prevents L and M inversion // +/- 10 channels prevents L and M inversion
fit->SetParLimits(1 + 2*i, fit->GetParameter(1+2*i)-10, fit->SetParLimits(1 + 2*i, fit->GetParameter(1+2*i)-dChannel,
fit->GetParameter(1+2*i)+10); fit->GetParameter(1+2*i)+dChannel);
fit->SetParLimits(2 + 2*i, 0, 1e5); fit->SetParLimits(2 + 2*i, 0, 1e5);
} }
// width of 975 keV line within 20% of 480 keV line // fix arbitrary amplitude for M component
// parameter is not used in the fit
fit->FixParameter(6, 50);
// width of 975 keV line constrainted wrt 480 keV line
double factor = 1.2; // 20%
if (p == 1) { if (p == 1) {
fit->SetParLimits(0, sigma480/1.1, sigma480*1.1); fit->SetParLimits(0, sigma480/factor, sigma480*factor);
} }
// fit spectrum // fit spectrum
...@@ -316,8 +336,14 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -316,8 +336,14 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
Double_t param[3]; Double_t param[3];
param[2] = fit->GetParameter(0); param[2] = fit->GetParameter(0);
for (UInt_t i = 0; i < peakListForFit.size(); ++i) { for (UInt_t i = 0; i < peakListForFit.size(); ++i) {
param[0] = fit->GetParameter(2 + 2*i); // param[0] = fit->GetParameter(2 + 2*i);
param[1] = fit->GetParameter(1 + 2*i); param[1] = fit->GetParameter(1 + 2*i);
if (i == 2) { // M component
param[0] = fit->GetParameter(2*i) / 4.;
}
else {
param[0] = fit->GetParameter(2 + 2*i);
}
signalFcn->SetParameters(param); signalFcn->SetParameters(param);
signalFcn->DrawCopy("same"); signalFcn->DrawCopy("same");
} }
...@@ -341,10 +367,10 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_ ...@@ -341,10 +367,10 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
/////////////////////////////////////////////////////////////////////////// ///////////////////////////////////////////////////////////////////////////
// fill fwhm data // fill fwhm data
if (p) { if (p) {
grFWHMhigh->SetPoint(n, n, fcalib->GetParameter(1)*2.35*fit->GetParameter(0)); grFWHMhigh->SetPoint(n, n, fcalibFull->GetParameter(1)*2.35*fit->GetParameter(0));
} }
else { else {
grFWHMlow ->SetPoint(n, n, fcalib->GetParameter(1)*2.35*fit->GetParameter(0)); grFWHMlow ->SetPoint(n, n, fcalibFull->GetParameter(1)*2.35*fit->GetParameter(0));
} }
} }
...@@ -449,7 +475,7 @@ void AddPeak(Double_t channel, Double_t cmin, Double_t cmax) ...@@ -449,7 +475,7 @@ void AddPeak(Double_t channel, Double_t cmin, Double_t cmax)
{ {
if (channel>cmin && channel<cmax) { if (channel>cmin && channel<cmax) {
peakListForFit.push_back(channel); peakListForFit.push_back(channel);
std::cout << "\t\t-> Added 1 peak to fit list at ch " << channel << "\n"; std::cout << "\t. Added 1 peak at ch " << channel << "\n";
} }
} }
...@@ -459,9 +485,6 @@ Double_t fpeaks(Double_t *x, Double_t *par) ...@@ -459,9 +485,6 @@ Double_t fpeaks(Double_t *x, Double_t *par)
{ {
Double_t result = 0; Double_t result = 0;
// linear background
// Double_t result = par[0] + par[1]*x[0];
// gaussian width // gaussian width
Double_t sigma = par[0]; Double_t sigma = par[0];
for (UInt_t p = 0; p < peakListForFit.size(); p++) { for (UInt_t p = 0; p < peakListForFit.size(); p++) {
...@@ -476,6 +499,37 @@ Double_t fpeaks(Double_t *x, Double_t *par) ...@@ -476,6 +499,37 @@ Double_t fpeaks(Double_t *x, Double_t *par)
Double_t fpeaks2(Double_t *x, Double_t *par)
{
Double_t result = 0;
// gaussian width
Double_t sigma = par[0];
// K component
result = par[2]*TMath::Gaus(x[0],par[1],sigma);
// L component
result += par[4]*TMath::Gaus(x[0],par[3],sigma);
// M component
result += par[4]/4*TMath::Gaus(x[0],par[5],sigma);
/*
for (UInt_t p = 0; p < peakListForFit.size(); p++) {
// case where brho is allowed to vary in fit
Double_t mean = par[2*p+1];
Double_t norm = par[2*p+2];
result += norm*TMath::Gaus(x[0],mean,sigma);
}
par[4] = 3.98*par[6];
*/
return result;
}
// Gaussian Peak function // Gaussian Peak function
Double_t gaussianPeak(Double_t *x, Double_t *par) Double_t gaussianPeak(Double_t *x, Double_t *par)
{ {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment