#include <iostream>
#include <string>
#include <sstream>
#include <fstream>
#include <vector>
#include "TH1.h"
#include "TF1.h"
#include "TSpectrum.h"
using namespace std;
string floatToStr(float f){
std::ostringstream os;
os << f ;
return (os.str());
}
// laboratory calibration
double labkevtochannel (double kev) {
double a = 0.000162834;
double b = 2.59414;
double c = -28.8893;
return -(b/(2*a))+sqrt(((pow(b,2))/(4*pow(a,2)))-(c-kev)/a);
}
double labchanneltokev (int channel) {
double a = 0.000162834;
double b = 2.59414;
double c = -28.8893;
88
APPENDIX B. PROGRAM CODE LISTING
return a*pow(channel,2)+b*channel+c;
}
double kevtochannel (double kev) {
return labkevtochannel(kev);
}
double channeltokev (int channel) {
return labchanneltokev(channel);
}
// Kista station calibration
double kistakevtochannel (double kev) {
double a = 0;
double b = 2.393936;
double c = -39.480137;
return (kev-c)/b;
}
double kistachanneltokev (int channel) {
double a = 0;
double b = 2.393936;
double c = -39.480137;
return b*channel+c;
}
double fwhmfromchannel (int channel) {
double a = 0.0921692;
double b = 20.6043;
double x = a*channel+b;
double y = channeltokev(channel);
return kevtochannel(y+(x/2.0)) - kevtochannel(y-(x/2.0));
}
double photoefficiency (int channel) {
// calculates how many of the
// emitted photons end up in the photo peak
double a = -0.000111023;
double b = 0.190384;
return a*channeltokev(channel)+b;
}
double photoefficiencyint (int channel) {
// calculates how many of the
// emitted photons end up in the photo peak
89
APPENDIX B. PROGRAM CODE LISTING
double a = -0.000111023;
double b = 0.190384;
return 0.5*a*channeltokev(channel)*channeltokev(channel)
+b*channeltokev(channel);
}
double comptonrelativeheight (double photopeakpos) {
// calculates the ratio between the integrated photopeak and
// the compton continuum height
// in case the Kista station
// data is used - conversion of the channels
photopeakpos = labchanneltokev(kevtochannel(photopeakpos));
double a = 1.16769e-005;
double b = 0.00302349;
double d = -1.12727e-008;
return d*pow(photopeakpos,3)
+a*pow(photopeakpos,2)
+b*photopeakpos;
}
double multiplecomptonrelativeheight (double photopeakpos) {
// calculates the ratio between
// the integrated photopeak and the height
// at the position between photo peak and compton edge
// in case the Kista station
// data is used - conversion of the channels
photopeakpos = labchanneltokev(kevtochannel(photopeakpos));
double a = 2.82596e-007;
double b = 0.000970274;
return a*photopeakpos+b;
}
double comptonedgerelativeheight (double photopeakpos) {
// calculates the ratio between the integrated photopeak
// and the height at the position of the compton edge
// in case the Kista station
// data is used - conversion of the channels
photopeakpos = labchanneltokev(kevtochannel(photopeakpos));
double a = 4.07404e-006;
double b = 0.00251556;
return a*photopeakpos+b;
}
90
APPENDIX B. PROGRAM CODE LISTING
double sigmafromchannel (int channel) {
return fwhmfromchannel(channel)/2.35;
}
int comptonedgepos(double peakpos) {
double Epeak = channeltokev(peakpos);
double tmp = 2*(Epeak/511);
double Eedge = Epeak*(tmp/(1.0+tmp));
return kevtochannel(Eedge);
}
int round(double x) {
if (floor(x) > 0.5) {
return ceil(x);
}
return floor(x);
}
string removePeak(double p,TH1I *c,
int minchannel,string outfile,bool backgroundappr = 0) {
cout << "Unfolding the peak at channel " << p << "\n";
// positions
double peakposexact = p;
int peakpos = round(peakposexact);
double backtmp = (2*channeltokev(peakposexact))/511.0;
double comptonedgepos =
kevtochannel(channeltokev(peakposexact)*(backtmp/(1+backtmp)));
cout << "Compton edge should be at channel "
<< comptonedgepos << "\n";
double backscatterpeakpos =
kevtochannel(channeltokev(peakpos)/(1+backtmp));
cout << "Backscatter peak should be at channel "
<< backscatterpeakpos << "\n";
// heights
int peakheight = c->GetBinContent(round(peakpos));
// search FWHM
int leftFWHM = peakpos;
int rightFWHM = peakpos;
int leftend,rightend;
double derthreshold = 0.05;
while (
c->GetBinContent(leftFWHM) >= (c->GetBinContent(peakpos)/2)
91
APPENDIX B. PROGRAM CODE LISTING
&& c->GetBinContent(rightFWHM) >= (c->GetBinContent(peakpos)/2)
&& leftFWHM > 0
&& rightFWHM < 1024) {
leftFWHM--;
rightFWHM++;
}
double sigma = (rightFWHM - leftFWHM - 1)/2.35;
cout << "Estimated FWHM of photopeak: "
<< (rightFWHM - leftFWHM - 1) << "\n";
cout << "Estimated sigma of photopeak: " << sigma << "\n";
// go to left
double maxder =
(((double) (c->GetBinContent(leftFWHM+2)
-c->GetBinContent(leftFWHM-2)))
/((double) 4));
for (int i = peakpos;
(i > leftFWHM
|| (((double) (c->GetBinContent(i+2)-c->GetBinContent(i-2)))
/((double) 4)) > (derthreshold*maxder))
&& i >= 0
&& i <= 1024;i--) {
double tmpder =
(((double) (c->GetBinContent(i+2)-c->GetBinContent(i-2)))
/((double) 4));
if (tmpder > maxder) {
maxder = tmpder;
}
}
leftend = i;
// go to right
double minder = (
((double) (c->GetBinContent(rightFWHM+2)
-c->GetBinContent(rightFWHM-2)))
/((double) 4));
for (int i = peakpos; (i < rightFWHM ||
(
((double) (c->GetBinContent(i+2)-c->GetBinContent(i-2)))
/((double) 4)) < (derthreshold*minder))
&& i >= 0
&& i <= 1024 ;i++) {
double tmpder = (
(
(double) (c->GetBinContent(i+2)-c->GetBinContent(i-2))
)
92
APPENDIX B. PROGRAM CODE LISTING
/((double) 4));
if (tmpder < minder) {
minder = tmpder;
}
}
rightend = i;
cout << "Estimated fitting area for peak: channel "
<< leftend << " to " << rightend << "\n";
// now the background estimate:
// a line through the left and right end of the peak
// f(x) = a*x + b
// calculate a and b
// (only if background approximation is applied)
double baka = 0;
double bakb = 0;
if (backgroundappr) {
baka = ((double) c->GetBinContent(rightend)
- c->GetBinContent(leftend))
/(rightend-leftend);
bakb = (double) c->GetBinContent(leftend) - baka*leftend;
cout << "Background function: f(x) = "
<< baka << "*x + "
<< bakb << "\n";
// create a double Gaussian for the left and right side
TF1 *faL = new TF1("faL","gaus(0)+[3]*x+[4]",0,1024);
TF1 *faR = new TF1("faR","gaus(0)+[3]*x+[4]",0,1024);
// set parameters for the background
faL->SetParameter(3,baka);
faL->SetParameter(4,bakb);
faR->SetParameter(3,baka);
faR->SetParameter(4,bakb);
faL->SetParLimits(3,baka,baka);
faL->SetParLimits(4,bakb,bakb);
faR->SetParLimits(3,baka,baka);
faR->SetParLimits(4,bakb,bakb);
} else {
cout << "No background approximation function applied!\n";
// create a double Gaussian for the left and right side
TF1 *faL = new TF1("faL","gaus(0)",0,1024);
TF1 *faR = new TF1("faR","gaus(0)",0,1024);
}
// Gaussian fit for left of the peak
faL->SetParLimits(0,peakheight/2.0,(3.0/2.0)*peakheight);
93
APPENDIX B. PROGRAM CODE LISTING
faL->SetParLimits(1,peakpos-1,peakpos+1);
faL->SetParLimits(2,sigma/2.0,(3.0/2.0)*sigma);
cout << "Fitting left half of the peak\n:";
c->Fit("faL","","",leftend,peakpos);
// Gaussian fit for right of the peak
faR->SetParLimits(0,peakheight/2.0,(3.0/2.0)*peakheight);
faR->SetParLimits(1,peakpos-1,peakpos+1);
faR->SetParLimits(2,sigma/2.0,(3.0/2.0)*sigma);
cout << "Fitting right half of the peak\n:";
c->Fit("faR","","",peakpos,rightend);
// after the fit both functions have to be
// reset so that they won’t use the background
// this is simply done by setting both
// parameters of the background line to zero
if (backgroundappr) {
faL->SetParameter(3,0);
faL->SetParameter(4,0);
faR->SetParameter(3,0);
faR->SetParameter(4,0);
}
// Furthermore we can readjust the boundaries of the peak
for (int i = peakpos; faR->Eval(i,0,0,0) > 1;i++) {
rightend = i;
}
for (int i = peakpos; faL->Eval(i,0,0,0) > 1;i--) {
leftend = i;
}
cout << "New estimate for peak width: channel "
<< leftend << " to " << rightend << "\n";
// giving an estimate for the height of the compton continuum
// first we integrate to see how many counts
// are in the photo peak
double photosum = 0;
for (int k = comptonedgepos;k <= 1024;k++) {
if (k <= peakpos) {
photosum += faL->Eval(k,0,0,0);
} else {
photosum += faR->Eval(k,0,0,0);
}
}
94
APPENDIX B. PROGRAM CODE LISTING
cout << "Integrated photo peak: " << photosum << "\n";
int comptonheight = round(
(photosum*comptonrelativeheight(peakposexact))
/comptonedgepos);
// the backscatter peak becomes
// much smaller with higher energies
double bskev1 = 661.62; // at the Cs-137 peak
double bsheight1 = 1.75;
// 175% of the compton continuum at Cs-137
double bskev2 = 1836; // at the Y-88 peak
double bsheight2 = 1.1;
// 110% of the compton continuum (estimate)
double bstmpa = (bsheight2-bsheight1)/(bskev1-bskev2);
double bstmpb = bsheight1 - bstmpa*bskev1;
int bspeakheight =
round(
(bstmpa*channeltokev(peakposexact)+bstmpb)
*comptonheight);
cout << "Estimated height of the compton continuum: "
<< comptonheight << "\n";
cout << "Estimated height of the backscatterpeak: "
<< bspeakheight << "\n";
int comptonedgelevel =
round(comptonedgerelativeheight(peakposexact)*photosum);
int multiplecomptonlevel =
round(multiplecomptonrelativeheight(peakposexact)*photosum);
// logistical function for the compton continuum
// 0: G
// 1: a
// 2: k
// 3: basic level multiple comptons
// between compton edge and photo peak
bool comptonactive = 1;
// indicates of the compton shift can be properly reduced
TF1 *continuum =
new TF1("continuum",
"[0]*(1-(1/(1+[1]*exp(-1.0*[2]*x))))+[3]",0,1024);
double G = comptonheight - multiplecomptonlevel;
if (comptonheight > 0) {
double b1 = ((double) comptonedgelevel - multiplecomptonlevel)
/comptonheight; // Coefficient 1: level at compton edge
double x1 = comptonedgepos;
double b2 = 1.0/G; // Coefficient 2:
95
APPENDIX B. PROGRAM CODE LISTING
// virtual zero-level between compton edge and peak position
double x2 = (x1+peakposexact)/2.0;
// equation for determination of (bx/(1-bx)) = a*exp(-k*bin)
cout << "G: " << G << "\n";
cout << "b1: " << b1 << "\n";
cout << "x1 (comptonedgepos): " << x1 << "\n";
cout << "b2: " << b2 << "\n";
cout << "x2 (between compton and peak): " << x2 << "\n";
if(G > 0 && b1 < 1 && b2 < 1 && b1 > 0 && b2 > 0 && x1 != x2) {
double contk = log((b1*(1-b2))/((1-b1)*b2))/(-x1+x2);
cout << "contk: " << contk << "\n";
double conta = (b1/(1-b1))*(1/exp(-contk*x1));
cout << "conta: " << conta << "\n";
continuum->SetParameter(0,G);
continuum->SetParameter(1,conta);
continuum->SetParameter(2,contk);
continuum->SetParameter(3,multiplecomptonlevel);
}
else {
comptonactive = 0; // compton shift can be properly reduced
}
} else {
comptonactive = 0; // compton shift can be properly reduced
}
//modification for elevated area
//of compton continuum near the compton edge
// we use a polynomial
// f(x) = a*x**2 + b*x + c with three characteristics
// channel zero (called k):
// elevated with half the elevation of the right half
// compton edge (called m):
// elevated with 10% at 661.62 keV and 70 % at 1836 keV
// middle of those two (called zero):
// no elevation, which leads to c = 0
double tmpq = 1.9595e-007*pow(channeltokev(peakposexact),2)
+ 2.14999e-005*channeltokev(peakposexact);
double tmpp = tmpq/2.0;
double tmpzero = comptonedgepos/2.0;
double tmpk = -tmpzero;
double tmpm = tmpzero;
double tmpa = (tmpq-(tmpp*tmpm)/tmpk)/(tmpm*tmpm-tmpk*tmpm);
double tmpb = (tmpp-tmpa*tmpk*tmpk)/tmpk;
96
APPENDIX B. PROGRAM CODE LISTING
TF1 *comptoncontmod1 =
new TF1("comptoncontmod1",
"1+[0]*(x-[2])**2+[1]*(x-[2])",0,1024);
comptoncontmod1->SetParameter(0,tmpa);
comptoncontmod1->SetParameter(1,tmpb);
comptoncontmod1->SetParameter(2,tmpzero);
// double Gaussian for the backscatter peak
bool backscatteractive = 1;
// indicates if backscatterpeak can properly be reduced
int backscatterpeakheight =
bspeakheight - (G+multiplecomptonlevel);
TF1 *backscatterL = new TF1("backscatterL","gaus(0)",0,1024);
TF1 *backscatterR = new TF1("backscatterR","gaus(0)",0,1024);
if (backscatterpeakheight > 1) {
double cL =
(-(1.0/4.0)*backscatterpeakpos)
/sqrt(-2*log(1.0/backscatterpeakheight));
double cR =
(kevtochannel(290)-backscatterpeakpos)
/sqrt(-2*log(1.0/backscatterpeakheight));
// parameters
// 0: backscatterpeakheight
backscatterL->SetParameter(0,backscatterpeakheight);
backscatterR->SetParameter(0,backscatterpeakheight);
// 1: backscatterpeakpos
backscatterL->SetParameter(1,backscatterpeakpos);
backscatterR->SetParameter(1,backscatterpeakpos);
// 2: cL/cR
backscatterL->SetParameter(2,cL);
backscatterR->SetParameter(2,cR);
} else {
backscatteractive = 0;
}
// escape functions
TF1 *singleescape = new TF1("backscatterR","gaus(0)",0,1024);
TF1 *doubleescape = new TF1("backscatterR","gaus(0)",0,1024);
singleescape->SetParameter(1,
kevtochannel(channeltokev(peakposexact)-511));
doubleescape->SetParameter(1,
kevtochannel(channeltokev(peakposexact)-2*511));
singleescape->SetParameter(2,
sigmafromchannel(kevtochannel(
channeltokev(peakposexact)-511)));
doubleescape->SetParameter(2,
97
APPENDIX B. PROGRAM CODE LISTING
sigmafromchannel(kevtochannel(
channeltokev(peakposexact)-2*511)));
if (channeltokev(peakposexact) > 2*511) {
cout << "Single Escape Peak:"
<< kevtochannel(channeltokev(peakposexact)-511) << "\n";
cout << "Double Escape Peak:"
<< kevtochannel(channeltokev(peakposexact)-2*511) << "\n";
double escapepeakheight =
(2450.0/20584839.0)
*(channeltokev(peakposexact)-2*511)*peakheight;
}
else {
double escapepeakheight = 0;
}
doubleescape->SetParameter(0,escapepeakheight);
singleescape->SetParameter(0,escapepeakheight);
// parts of the function
int gauss,logist,backscatterbin,singleescapebin,doubleescapebin;
// dampening estimate
int sumin,sumrem; // sums of initial and remaining values
// open output file
ofstream offile;
offile.open (outfile.c_str(), ios::out);
// run through and deduct
for (int i = minchannel; i <= 1024;i++) {
// save initial value
int initial = c->GetBinContent(i);
// set all parts to zero
gauss = logist = backscatterbin
= singleescapebin = doubleescapebin = 0;
// Gaussian peak
if (i < p) {
int gauss = round(faL->Eval(i,0,0,0));
} elseif (i == p) {
int gauss = round((faL->Eval(i,0,0,0)
+ faR->Eval(i,0,0,0))/2.0);
} else {
int gauss = round(faR->Eval(i,0,0,0));
}
98
APPENDIX B. PROGRAM CODE LISTING
// continuum
if (
(continuum->Eval(i,0,0,0)
* comptoncontmod1->Eval(i,0,0,0)
) >
faL->Eval(i,0,0,0)
&& i < p && comptonactive) {
int logist = round(continuum->Eval(i,0,0,0)
* comptoncontmod1->Eval(i,0,0,0));
int gauss = 0;
} else {
int logist = 0;
}
// backscatter peak
if (backscatteractive) {
if (i < backscatterpeakpos) {
int backscatterbin = round(backscatterL->Eval(i,0,0,0));
} elseif (i == backscatterpeakpos) {
int backscatterbin = round((backscatterL->Eval(i,0,0,0)
+ backscatterR->Eval(i,0,0,0))/2.0);
} else {
int backscatterbin = round(backscatterR->Eval(i,0,0,0));
}
} else {
int backscatterbin = 0;
}
// escape peaks
if(channeltokev(peakposexact)-2*511 > 0) {
int singleescapebin = round(singleescape->Eval(i,0,0,0));
int doubleescapebin = round(doubleescape->Eval(i,0,0,0));
}
// subtract the parts
c->SetBinContent(i,
initial-gauss-logist-
backscatterbin-singleescapebin-doubleescapebin);
// the protocol file
offile << i << "\t"; // 1: channel number
offile << channeltokev(i) << "\t"; // 2: corresponding energy
offile << initial << "\t"; // 3: initial value
offile << gauss << "\t"; // 4: Photo peak
offile << baka*i+bakb << "\t";
99
APPENDIX B. PROGRAM CODE LISTING
// 5: value of the background approximation
// function at this position
offile << logist << "\t"; // 6: Logistic function
offile << backscatterbin << "\t"; // 7: Backscatter peak
offile << singleescapebin << "\t"; // 8:Single Escape Peak
offile << doubleescapebin << "\t"; // 9: Double Escape Peak
offile
<< (gauss+logist+backscatterbin+
singleescapebin+doubleescapebin)
<< "\t";
// 10: the whole modelled response function combined
offile << c->GetBinContent(i) << "\t"; // 11: remainder
offile << endl;
// estimate for the dampening
sumin += abs(initial);
sumrem += abs(c->GetBinContent(i));
}
cout << "Initial sum:\t" << sumin << "\n";
cout << "Remaining sum:\t" << sumrem << "\n";
cout << "Dampening:\t"
<< 100-(sumrem/(double) sumin)*100 << "% \n";
string logstring;
logstring.append("Counts in photo peak: ");
logstring.append(floatToStr(photosum));
logstring.append("\n");
return logstring;
offile.close();
}
int test5c ()
{
gROOT->Reset();
int channels = 1024;
int minchannel = 27;
TH1I *counts =
new TH1I("Histogram of Counts","Counts",channels,0,channels);
counts->Reset();
100
APPENDIX B. PROGRAM CODE LISTING
// readfile
ifstream in;
cout << "Choose one of the samples." << endl;
cout << "Predefined peak positions:" << endl;
cout << "(1) Cs-137 (point)" << endl;
cout << "(2) Cs-137 (extended)" << endl;
cout << "(3) Co-60" << endl;
cout << "(4) Y-88" << endl;
cout << "(5) Ba-133" << endl;
cout << "(6) Na-22" << endl;
cout << "Automatic peak search:" << endl;
cout << "(11) Cs-137 (point)" << endl;
cout << "(12) Cs-137 (extended)" << endl;
cout << "(13) Co-60" << endl;
cout << "(14) Y-88" << endl;
cout << "(15) Ra-226" << endl;
cout << "(16) Ba-133" << endl;
cout <<
"(21) Kista station spectrum with extended Cs-137 source"
<< endl;
cout << "(22) Real spectrum (1 hour)" << endl;
cout << "(23) Real spectrum (1 hour, with Co-60 source)"
<< endl;
int a;
cin >> a;
double x,y;
string datfile;
switch (a) {
case 1: in.open("C:/root/macros/histocs137.dat"); break;
case 2: in.open("C:/root/macros/histocs137ex.dat"); break;
case 3: in.open("C:/root/macros/histoco60.dat"); break;
case 4: in.open("C:/root/macros/histoy88.dat"); break;
case 5: in.open("C:/root/macros/histoba133.dat"); break;
case 6: in.open("C:/root/macros/histona22.dat"); break;
case 11: in.open("C:/root/macros/histocs137.dat"); break;
case 12: in.open("C:/root/macros/histocs137ex.dat"); break;
case 13: in.open("C:/root/macros/histoco60.dat"); break;
case 14: in.open("C:/root/macros/histoy88.dat"); break;
case 15: in.open("C:/root/macros/histora226.dat"); break;
case 16: in.open("C:/root/macros/histoba133.dat"); break;
case 21: in.open("C:/root/macros/historoof1.dat"); break;
case 22: in.open("C:/root/macros/historealnormal.dat"); break;
case 23:
in.open("C:/root/macros/historealcontaminated.dat"); break;
101
APPENDIX B. PROGRAM CODE LISTING
}
if (!in.is_open()) { printf("error!"); }
int i = 1;
while(in.good()) {
in >> x >> y;
i++;
if (x < minchannel) {
y = 0;
}
counts->SetBinContent(x,y);
}
in.close();
int sum1;
for (int i = minchannel;i <= 1024;i++) {
sum1 += abs(counts->GetBinContent(i));
}
TSpectrum *spec = new TSpectrum(1024);
switch (a) {
case 1:
removePeak(262.75,counts,minchannel,
"C:/root/macros/outputcs137.dat");
break;
case 2:
removePeak(262.75,counts,minchannel,
"C:/root/macros/outputcs137ex.dat");
break;
case 3:
removePeak(510.0,counts,minchannel,
"C:/root/macros/outputco60-510.dat");
removePeak(450.0,counts,minchannel,
"C:/root/macros/outputco60-450.dat");
break;
case 4:
removePeak(689.0,counts,minchannel,
"C:/root/macros/outputy88-689.dat");
removePeak(351.0,counts,minchannel,
"C:/root/macros/outputy88-351.dat");
break;
case 5:
removePeak(kevtochannel(383.851),counts,minchannel,
"C:/root/macros/outputba133-384.dat");
removePeak(kevtochannel(356.017),counts,minchannel,
"C:/root/macros/outputba133-356.dat");
102
APPENDIX B. PROGRAM CODE LISTING
removePeak(kevtochannel(302.853),counts,minchannel,
"C:/root/macros/outputba133-302.dat");
removePeak(kevtochannel(276.398),counts,minchannel,
"C:/root/macros/outputba133-276.dat");
break;
case 6:
removePeak(kevtochannel(1274.53),counts,minchannel,
"C:/root/macros/outputna22-1275.dat");
removePeak(kevtochannel(511),counts,minchannel,
"C:/root/macros/outputna22-511.dat");
break;
case 7:
removePeak(kevtochannel(383.851),counts,minchannel,
"C:/root/macros/outputba133-384.dat");
removePeak(kevtochannel(356.017),counts,minchannel,
"C:/root/macros/outputba133-356.dat");
removePeak(kevtochannel(302.853),counts,minchannel,
"C:/root/macros/outputba133-302.dat");
removePeak(kevtochannel(276.398),counts,minchannel,
"C:/root/macros/outputba133-276.dat");
break;
case 8:
removePeak(kevtochannel(383.851),counts,minchannel,
"C:/root/macros/outputba133-384.dat");
removePeak(kevtochannel(356.017),counts,minchannel,
"C:/root/macros/outputba133-356.dat");
removePeak(kevtochannel(302.853),counts,minchannel,
"C:/root/macros/outputba133-302.dat");
removePeak(kevtochannel(276.398),counts,minchannel,
"C:/root/macros/outputba133-276.dat");
break;
default:
std::ostringstream logstream;
vector<double> v;
logstream << "Report about dampening" << endl;
int xpeaktopheight;
float xpeaktoppos;
bool flag = 1;
// overload of calibration functions for Kista station samples
if (a > 20) {
int kevtochannel (double kev) {
return kistakevtochannel(kev);
}
double channeltokev (int channel) {
return kistachanneltokev(channel);
103
APPENDIX B. PROGRAM CODE LISTING
}
}
for (int l = 1;l <= 2 && flag;l++) {
//search peaks
Int_t nfound = spec->Search(counts,2);
cout << "Found " << nfound << " candidate peaks to fit\n";
Float_t *xpeaks = spec->GetPositionX();
bool xpeakfound = 0;
xpeaktopheight = 0;
float xpeaktmp;
// sort the peaks by position
while (flag) {
flag = 0;
for (int i = 1; i < nfound; i++) {
if (xpeaks[i-1] < xpeaks[i]) {
xpeaktmp = xpeaks[i];
xpeaks[i] = xpeaks[i-1];
xpeaks[i-1] = xpeaktmp;
flag = 1;
}
}
}
// list the found peaks
for (int i = 0; i < nfound; i++) {
cout << "Channel " << xpeaks[i]
<< " (" << channeltokev(xpeaks[i])
<< " keV)" << endl;
}
flag = 0;
for (int i = 0; i < nfound; i++) {
// check if this peak has already been reduced
vector<double>::iterator iter1 = v.begin();
bool doubleredflag = 0;
while( iter1 != v.end()) {
if (*iter1 == xpeaks[i]) {
doubleredflag = 1;
}
iter1++;
}
if (counts->GetBinContent(round(xpeaks[i])) > 0
&& channeltokev(xpeaks[i]) > 100
&& !doubleredflag) {
// run peak dampening
104
APPENDIX B. PROGRAM CODE LISTING
logstream << "Reduction round: " << l << endl;
logstream << "Dampening at channel "
<< xpeaks[i] << " ("
<< channeltokev(xpeaks[i])
<< " keV)" << endl;
string basestr = "C:/root/macros/outputautomatic-";
basestr.append(floatToStr(xpeaks[i]));
basestr.append(".dat");
logstream << removePeak(xpeaks[i],counts,minchannel,basestr,1);
v.push_back(xpeaks[i]);
logstream << "----\n";
flag = 1;
}
}
}
cout << logstream.str();
break;
}
int sum2;
for (int i = minchannel;i <= 1024;i++) {
sum2 += abs(counts->GetBinContent(i));
}
if (sum1 == 0) {
cout << "no dampening possible" << endl;
} else {
cout << "Final dampening: "
<< (100-(sum2/(float) sum1)*100)
<< "% \n";
}
counts->Draw();
}