diff --git a/Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C b/Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C
index 0cd2ac7029635f51b9bae0658ecdac858d3f55de..688c770f47d4134f997959affef24c89793417dd 100644
--- a/Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C
+++ b/Projects/ComptonTelescope/online/Calibration_DSSSD/Analyse207Bi.C
@@ -31,6 +31,7 @@
#include "TLine.h"
#include "TError.h"
#include "TROOT.h"
+#include "TMath.h"
// C++ headers
#include <iostream>
@@ -54,6 +55,7 @@ static std::vector<Double_t> peakListForFit;
// functions
void AddPeak(Double_t, 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*);
@@ -89,7 +91,7 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
///////////////////////////////////////////////////////////////////////////
// open output pdf file
string label = sname.substr(4, 6);
- label += sname.substr(16, 20);
+ label += sname.substr(16, 21);
label += (isPside) ? "_pside" : "_nside";
TString pdfname = Form("Analyse207Bi_%s.pdf", label.c_str());
can->Print(Form("%s[", pdfname.Data()));
@@ -101,7 +103,8 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
///////////////////////////////////////////////////////////////////////////
// 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
@@ -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
// treat all channels or specified one
if (channel < 0 || channel == n) {
- cout << "Analysing channel " << n << "\n";
+ cout << "\n--------------- Analysing channel " << n << " ---------------\n";
///////////////////////////////////////////////////////////////////////////
// get histogram
@@ -168,8 +171,8 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
auto s = new TSpectrum();
Int_t npeaks = s->Search(h, 6, "", 0.6);
h->DrawCopy();
- if (npeaks == 1) cout << "pedestal found!\n";
- else cout << "PROBLEM with pedestal!\n";
+ if (npeaks == 1) cout << ". Pedestal found -> OK\n";
+ else cout << ". Pedestal PROBLEM\n";
Double_t *xpeaks = s->GetPositionX();
// fit
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_
// search peaks after background subtraction
npeaks = s->Search(h, 6, "", 0.3);
h->DrawCopy();
- cout << "\t-> found " << npeaks << " peaks";
+ cout << ". " << npeaks << " main transitions found ";
if (npeaks == 2) cout << "\t-> OK!\n";
else cout << "\t-> PROBLEM!\n";
// get peak position
@@ -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
// uses calibration parameters from previous strip calibration
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;
if (energy > mainLineEnergy[0]+200) {
- energyGuess = (mainLineEnergy[0]-fcalib->GetParameter(0))/fcalib->GetParameter(1);;
+ energyGuess = (mainLineEnergy[0]-fcalibFull->GetParameter(0))/fcalibFull->GetParameter(1);;
}
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);
}
@@ -230,9 +233,9 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
///////////////////////////////////////////////////////////////////////////
// rough energy calibration
auto grcalib = new TGraph(peakList.size(), &peakList[0], &mainLineEnergy[0]);
- grcalib->Fit("fcalib", "Q0");
- grOffset->SetPoint(n, n, fcalib->GetParameter(0));
- grGain ->SetPoint(n, n, fcalib->GetParameter(1)*10);
+ grcalib->Fit("fcalibFull", "Q0");
+ grOffset->SetPoint(n, n, fcalibFull->GetParameter(0));
+ grGain ->SetPoint(n, n, fcalibFull->GetParameter(1)*10);
// loop on "480"- and "975"-keV features
// p = 0 -> 480 keV
@@ -242,15 +245,22 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
can->cd(3+p);
pad = (TPad*) can->FindObject(Form("can_%d", 3+p));
// determine range for fitting based on rough calibration
- Double_t cmin = (mainLineEnergy[p]-40 -fcalib->GetParameter(0)) / fcalib->GetParameter(1);
- Double_t cmax = (mainLineEnergy[p]+120-fcalib->GetParameter(0)) / fcalib->GetParameter(1);
- std::cout << cmin << "\t" << cmax << "\n";
+ Double_t cmin = (mainLineEnergy[p]-40 -fcalibFull->GetParameter(0)) / fcalibFull->GetParameter(1);
+ Double_t cmax = (mainLineEnergy[p]+130-fcalibFull->GetParameter(0)) / fcalibFull->GetParameter(1);
+ 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->DrawCopy();
// search satellite peaks
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";
else cout << "\t-> PROBLEM!\n";
// get peak position
@@ -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
// should be removed
if (npeaks == 3) peakListForFit.erase(peakListForFit.begin()+1);
- // add CE-M component
- Double_t channel = (MLineEnergy[p] -fcalib->GetParameter(0)) / fcalib->GetParameter(1);
+ // add CE-M component based on local linear fit
+ 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);
// define parameter list for fit
std::vector<Double_t> paramList;
// width (sigma); common for all states
- paramList.push_back(4);
+ paramList.push_back(6);
// position and amplitude
for (UInt_t i = 0; i < peakListForFit.size(); ++i) {
paramList.push_back(peakListForFit[i]);
@@ -284,20 +298,26 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
}
// 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(¶mList[0]);
// positive amplitudes and position in range
+ int dChannel = 5;
for (UInt_t i = 0; i < peakListForFit.size(); ++i) {
// +/- 10 channels prevents L and M inversion
- fit->SetParLimits(1 + 2*i, fit->GetParameter(1+2*i)-10,
- fit->GetParameter(1+2*i)+10);
+ fit->SetParLimits(1 + 2*i, fit->GetParameter(1+2*i)-dChannel,
+ fit->GetParameter(1+2*i)+dChannel);
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) {
- fit->SetParLimits(0, sigma480/1.1, sigma480*1.1);
+ fit->SetParLimits(0, sigma480/factor, sigma480*factor);
}
// fit spectrum
@@ -316,8 +336,14 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
Double_t param[3];
param[2] = fit->GetParameter(0);
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);
+ if (i == 2) { // M component
+ param[0] = fit->GetParameter(2*i) / 4.;
+ }
+ else {
+ param[0] = fit->GetParameter(2 + 2*i);
+ }
signalFcn->SetParameters(param);
signalFcn->DrawCopy("same");
}
@@ -341,10 +367,10 @@ void Analyse207Bi(const char* name = "bb7_3309-7_bi207_20210126_13h09_run5_conv_
///////////////////////////////////////////////////////////////////////////
// fill fwhm data
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 {
- 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)
{
if (channel>cmin && channel<cmax) {
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)
{
Double_t result = 0;
- // linear background
- // Double_t result = par[0] + par[1]*x[0];
-
// gaussian width
Double_t sigma = par[0];
for (UInt_t p = 0; p < peakListForFit.size(); p++) {
@@ -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
Double_t gaussianPeak(Double_t *x, Double_t *par)
{