/*
XECON - a tool for analysis of radioxenon data
Copyright (C) 2009-2010 Anders Ringbom, Swedish Defence Research Agency (FOI)
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see .
*/
#include
//
//Define global constants
//
double ks = 1.6449;
double Xe133_lam=1.52956e-6; //Decay constant for 133Xe(g) (1/s)
double Xe131m_lam=6.72467e-7; //Decay constant for 131Xe(m)
double Xe133m_lam=3.66326e-6; //Decay constant for 133Xe(m)
double Xe135_lam=2.10657e-5; //Decay constant for 135Xe(g)
double bg_br_133Xe_80=0.37289; //Branching ratio for beta/80 keV gamma decay of 133Xe
double bg_br_133Xe_30=0.48940; //Branching ratio for beta/30 keV X-ray decay of 133Xe
double bg_br_131mXe_30=0.54000; //Branching ratio for EC/30 keV X-ray decay of 131mXe
double bg_br_133mXe_30=0.55770; //Branching ratio for EC/30 keV X-ray decay of 133mXe
double bg_br_135Xe=0.90000; //Branching ratio for beta/250 keV gamma decay of 135Xe
double Xe_air=0.087; //Volume (ml) of Xe per m3 of air
double Rn222_lam = 2.0974e-6; // radon decay constant
double bg_br_222Rn = 0.376; // radon branching ratio for 352 keV (ROI1)
double pi = 3.141592654;
#define MAXP 5000
#define MAXY 2100
int datatype(char *file, char *dtype){
char line[2000];
char br[100];
FILE *fp;
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"DATA_TYPE") == 0) {
//printf("%s\n", br);
sscanf(&line[0],"%s %s", &br, dtype);
}
}
fclose(fp);
return 0;
}
int header(char *file, char *station, char *detector){
char line[200];
char br[100];
FILE *fp;
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
//line[0]="";
while (fgets(line,100,fp)){
sscanf(line,"%s", br);
if ( strcmp(br,"#Header") == 0) {
//printf("%s\n", br);
fgets(line,100,fp);
sscanf(line,"%s %s",station,detector);
}
}
fclose(fp);
return 0;
}
double collection(char *file, char *start_date, char *start_time, char *stop_date, char *stop_time){
char line[2000];
char br[100];
FILE *fp;
double airvol;
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
double failurereturn=-1;
return failurereturn;
}
// line[0]=' ';
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#Collection") == 0) {
//printf("%s\n", br);
fgets(line,100,fp);
sscanf(line,"%s %s %s %s %lf",start_date, start_time, stop_date, stop_time, &airvol);
}
}
fclose(fp);
return airvol;
}
int processing(char *file, double* vol, double* yield){
char line[2000];
char br[100];
FILE *fp;
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#Processing") == 0) {
//printf("%s\n", br);
fgets(line,100,fp);
sscanf(line,"%lf %lf",vol, vol+1);
fgets(line,100,fp);
sscanf(line,"%lf %lf",yield, yield+1);
}
}
fclose(fp);
return 0;
}
int acquisition(char *file, char *start_date, char *start_time, double* atime){
char line[2000];
char br[100];
FILE *fp;
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
line[0]=' ';
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#Acquisition") == 0) {
//printf("%s\n", br);
fgets(line,100,fp);
sscanf(line,"%s %s %lf %lf",start_date, start_time, atime, atime+1);
}
}
fclose(fp);
return 0;
}
int g_energy(char *file, double* ene, double* chan, double* err){
char line[2000];
char br[100];
FILE *fp;
int n,np;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
n=0;
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
//printf("%s\n", br);
if ( strcmp(br,"#g_Energy") == 0) {
//printf("%s\n", br);
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%lf %lf %lf", ene+n,chan+n,err+n);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n++;
}
}
}
fclose(fp);
//printf("n gamma: %d\n",n);
return n;
}
int b_energy(char *file, double* ene, double* chan, double* err){
char line[2000];
char br[100],dummy[10];
FILE *fp;
int n,np;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#b_Energy") == 0) {
//printf("%s\n", br);
n = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%lf %s %lf %lf", ene+n,dummy,chan+n,err+n);
//printf("%s\n",line);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n++;
}
}
}
fclose(fp);
//printf("n beta: %d\n",n);
return n;
}
int g_resolution(char *file, double* ene, double* res, double* err){
char line[2000];
char br[100];
FILE *fp;
int n,np;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#g_Resolution") == 0) {
//printf("%s\n", br);
n = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%lf %lf %lf", ene+n,res+n,err+n);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n++;
}
}
}
fclose(fp);
return n;
}
int b_resolution(char *file, double* ene, double* res, double* err){
char line[2000];
char br[100];
FILE *fp;
int n,np;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#b_Resolution") == 0) {
//printf("%s\n", br);
n = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%lf %lf %lf", ene+n,res+n,err+n);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n++;
}
}
}
fclose(fp);
return n;
}
int ROI_limits(char *file, int* roi, double* b_low, double* b_high, double* g_low, double* g_high){
char line[2000];
char br[100];
FILE *fp;
int n,np,m;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#ROI_Limits") == 0) {
// printf("%s\n", br);
n = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%d %lf %lf %lf %lf",roi+n, b_low+n,b_high+n,g_low+n,g_high+n);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n++;
}
}
}
fclose(fp);
return n;
}
int GetGate(char *file, double* coord){
char line[2000];
char br[100];
FILE *fp;
int n,np;
int m=0;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#Gate") == 0) {
//printf("%s\n", br);
n = 0;
m = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%lf %lf",coord+n, coord+n+1);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n=n+2;
m++;
}
}
}
fclose(fp);
return m;
}
void GetTweak(char *file, double* shift){
char line[2000];
char br[100];
FILE *fp;
int n,np;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#Tweak") == 0) {
//printf("%s\n", br);
n = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%lf",shift+n);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n=n+1;
}
}
}
fclose(fp);
}
double GetMemf(char *file){
char line[2000];
char br[100];
double mem;
FILE *fp;
int n,np;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#Memf") == 0) {
//printf("%s\n", br);
n = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%lf",&mem);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n=n+1;
}
}
}
fclose(fp);
return mem;
}
void GetFilenames(char* file, char* sf, char* gf, char* df){
char sfile1[100],gfile1[100],dfile1[100];
char datadir[150], fl[5];
char iline[400];
FILE *finput = fopen(file,"r");
fgets(iline,100,finput);
sscanf(iline,"%s",datadir);
while (fgets(&iline[0],300,finput)){
sscanf(iline,"%s\t%s\t%s\t%s",sfile1,gfile1,dfile1,fl);
if ( strcmp(fl,"x") == 0){
sprintf(sf,"%s%s",datadir,sfile1);
sprintf(gf,"%s%s",datadir,gfile1);
sprintf(df,"%s%s",datadir,dfile1);
strcpy(fl,"");
}
}
fclose(finput);
//printf(" leaving getfilenames\n");
}
int bg_efficiency(char *file, int* roi, double* eff, double* err){
char line[2000];
char br[100],dum[30];
FILE *fp;
int n,np,m;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#b-gEfficiency") == 0) {
//printf("%s\n", br);
n = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%s %d %lf %lf",dum,roi+n,eff+n,err+n);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n++;
}
}
}
fclose(fp);
return n;
}
int ratios(char *file, int* roi1, int* roi2, double* rat, double* rat_err){
char line[2000];
char br[100],dum[30];
FILE *fp;
int n,np,m;
char nu = ' ';
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return failurereturn;
}
while (fgets(line,100,fp)){
sscanf(line,"%s", &br);
if ( strcmp(br,"#Ratios") == 0) {
//printf("%s\n", br);
n = 0;
fgets(line,100,fp);
while ( nu != '#' ){
sscanf(line,"%s %d %d %lf %lf",dum,roi1+n, roi2+n,rat+n,rat_err+n);
fgets(line,100,fp);
sscanf(line,"%c", &nu);
n++;
}
}
}
fclose(fp);
return n;
}
int* g_spectrum(char *file, int* chan){
char line[2000];
char br[100],dum[100];
FILE *fp;
int n,np,m;
char nu = ' ';
int *hist;
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return &failurereturn;
}
while (fgets(line,100,fp)){
sscanf(&line[0],"%s", &br);
if ( strcmp(br,"#g_Spectrum") == 0) {
//printf("%s\n", br);
fgets(line,100,fp);
sscanf(line,"%d %d", chan, chan+1);
n = 0;
hist=(int*)calloc(chan[0],sizeof(int));
while (n < chan[0]) {
fgets(line,100,fp);
sscanf(line,"%d %d %d %d %d %d ", &m, hist+n, hist+n+1, hist+n+2, hist+n+3, hist+n+4);
n = n+ 5;
}
fgets(line,100,fp);
sscanf(line,"%c", &nu);
}
}
fclose(fp);
return hist;
}
int* b_spectrum(char *file, int* chan){
char line[2000];
char br[100],dum[100];
FILE *fp;
int n,np,m;
char nu = ' ';
int *hist;
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return &failurereturn;
}
while (fgets(line,100,fp)){
sscanf(&line[0],"%s", &br);
if ( strcmp(br,"#b_Spectrum") == 0) {
//printf("%s\n", br);
fgets(line,100,fp);
sscanf(line,"%d %d", chan, chan+1);
n = 0;
hist=(int*)calloc(chan[0],sizeof(int));
while (n < chan[0]) {
fgets(line,100,fp);
sscanf(line,"%d %d %d %d %d %d ", &m, hist+n, hist+n+1, hist+n+2, hist+n+3, hist+n+4);
n = n+ 5;
}
fgets(line,100,fp);
sscanf(line,"%c", &nu);
}
}
fclose(fp);
return hist;
}
int* bg_spectrum(char *file, int* chan){
char line[2000];
char br[100],dum[100];
FILE *fp;
int n,np,m,lpos;
char nu = ' ';
int *hist;
if(!(fp = fopen(file,"r"))){
printf("\nfailed to open\n%s\n",file);
int failurereturn=-1;
return &failurereturn;
}
while (fgets(line,100,fp)){
sscanf(&line[0],"%s", &br);
if ( strcmp(br,"#Histogram") == 0) {
//printf("%s\n", br);
fgets(line,100,fp);
sscanf(line,"%d %d %d %d", chan, chan+1, chan+2, chan+3);
n = 0;
m = 0;
hist=(int*)calloc(chan[0]*chan[1],sizeof(int));
char snm[20];
while (n < chan[0] ) {
fgets(line,2000,fp);
lpos = 0;
for(m=0;m 0; count--,i++)
{
f >> hist[i];
}
f.close();
return hist;
}
}
}
// Creates a projection of a 2D histogram using given limits in channels
//
void project(char* hname ,int b_low, int b_high, int g_low, int g_high){
char textg[50], textb[50];
gDirectory->Delete("hbg_projection_x");
gDirectory->Delete("hbg_projection_y");
TH2F *hist=(TH2F*)gDirectory->Get(hname);
hist->ProjectionY("gamma_proj", b_low, b_high);
hist->ProjectionX("beta_proj", g_low, g_high);
TCanvas *My_projections = new TCanvas("Projections",hname,600,600);
My_projections->Divide(1,2);
My_projections->cd(1);
gamma_proj->Draw();
TLatex *text = new TLatex();
text->SetNDC();
text -> SetTextSize(0.07);
sprintf(textg,"Beta limits (channels): %d - %d",b_low,b_high);
text -> DrawLatex(0.2, 0.92, textg);
My_projections->cd(2);
beta_proj->Draw();
sprintf(textb,"Gamma limits (channels): %d - %d",g_low,g_high);
text -> DrawLatex(0.2, 0.92, textb);
}
// Fills a 1D histogram with pulse heigh data from an IMS file
// The histogram should be defined as follows:
// TH1F *hist = new TH1F("hist","", ch_en[0]+1, 0, ch_en[0]+1);
//
void fill_1D_hist(int* val, int* ch_en, char *hname){
int x,N;
TH1F *hist=(TH1F*)gDirectory->Get(hname);
// Fill histogram
for (x=1 ;x < ch_en[0]; x++){
if(N=val[x]){
hist->Fill(x,N);
}
}
}
// Fills a 2D histogram with pulse heigh data
// The histogram should be defined as follows:
// TH2F *hist = new TH2F("hname","", ch_en[1]+1, 0, ch_en[1]+1,ch_en[0]+1, 0, ch_en[0]+1);
//
void fill2Dhist(int* val,int* ch_en, char *hname){
int x,y,N;
TH2F *hist=(TH2F*)gDirectory->Get(hname);
// Fill histogram
for (y=1 ; y < ch_en[0]; y++) { // loop over rows (gamma)
for (x=1 ;x < ch_en[1]; x++){ // loop over elements in row (beta)
if(N=val[ch_en[1]*y+x]){
hist->Fill(x,y,N);
}
}
}
}
// Fills a 2D histogram with pulse heigh data
// The histogram should be defined as follows:
// TH2F *hist = new TH2F("hname","", ch_en[1]+1, 0, ch_en[1]+1,ch_en[0]+1, 0, ch_en[0]+1);
//
void fill2Dhist_all_chan(int* val,int* ch_en, char *hname){
int x,y,N;
TH2F *hist=(TH2F*)gDirectory->Get(hname);
// Fill histogram
for (y=0 ; y < ch_en[0]; y++) { // loop over rows (gamma)
for (x=0 ;x < ch_en[1]; x++){ // loop over elements in row (beta)
if(N=val[ch_en[1]*y+x]){
hist->Fill(x,y,N);
}
}
}
}
// Fills a 2D histogram with pulse heigh data using a cut
// The histogram should be defined as follows:
// TH2F *hist = new TH2F("hname","", ch_en[1]+1, 0, ch_en[1]+1,ch_en[0]+1, 0, ch_en[0]+1);
//
void fill_cut_2Dhist(int* val,int* ch_en, char *hname, TCutG *gate){
int x,y,N;
TH2F *hist=(TH2F*)gDirectory->Get(hname);
// Fill histogram
for (y=0 ; y < ch_en[0]; y++) { // loop over rows (gamma)
for (x=0 ;x < ch_en[1]; x++){ // loop over elements in row (beta)
if(N=val[ch_en[1]*y+x]){
if(gate->IsInside(x,y+1)) hist->Fill(x,y,N);
}
}
}
}
// Adds histogram array val2 to array val1:
//
void histarr_add(int* val1, int* val2, int* ch_en){
int x,y,N;
for (y=0 ; y < ch_en[0]; y++) { // loop over rows (gamma)
for (x=0 ;x < ch_en[1]; x++){ // loop over elements in row (beta)
if(N=val2[ch_en[1]*y+x]){
val1[ch_en[1]*y+x] = val1[ch_en[1]*y+x] + val2[ch_en[1]*y+x];
}
}
}
}
//Adds 2D beta-gamma data from to IMS files and writes the result to a third file as a #Histogram block.
//
//Parameters:
//
//char *file1: name of file 1
//char *file2: name of file 2
//char *file1: name of file to store the resulting #Histogram block
//
void add_write_bg(char *file1, char *file2, char *sumfile){
//char *dfile1="d2_sum.phd";
//char *dfile2="d2_3.phd";
//char *sum_dfile="d2_sum.phd";
//
//Create and fill 2D histograms
//
int chan_ene_1[4], *bg_hist_1;
int chan_ene_2[4], *bg_hist_2;
bg_hist_1=bg_spectrum(file1, chan_ene_1);
bg_hist_2=bg_spectrum(file2, chan_ene_2);
histarr_add(bg_hist_1,bg_hist_2,chan_ene_1);
FILE *fout= fopen(sumfile,"w");
fprintf(fout,"#Histogram \n");
fprintf(fout,"%d %d %d %d\n", chan_ene_1[0],chan_ene_1[1],chan_ene_1[0],chan_ene_1[1]);
int x,y;
for (y=0 ; y < chan_ene_1[0]; y++) { // loop over rows (gamma)
for (x=0 ;x < chan_ene_1[1]; x++){ // loop over elements in row (beta)
fprintf(fout,"%d ",bg_hist_1[chan_ene_1[1]*y+x]);
}
fprintf(fout,"\n");
}
fclose(fout);
}
//
// returns bin with max value in range lowx - highx
//
int GetMaxBinRange(char *hname, int lowx, int highx){
TH2F *hist=(TH2F*)gDirectory->Get(hname);
int nb = hist->GetXaxis()->GetNbins();
int max, x, maxb;
max = 0; maxb = 0;
for (i=lowx;i<=highx;i++){
x = hist->GetBinContent(i);
if(x > max) {
max = x;
maxb = i;
}
}
return maxb;
}
void MakeGate(int np, double *gate, TCutG *mgate[30]){
int m=0; int gate1,gate2;
mgate[1] = new TCutG("ugate",np+1);
for (int i = 0; i <= 2*(np-1); i=i+2){
mgate[1]->SetPoint(m,gate[i],gate[i+1]);
// printf("m %d gate %f %f\n",m,gate[i],gate[i+1]);
m++;
}
mgate[1]->SetPoint(np,gate[0],gate[1]);
//printf("np %d\n",np);
}
// Performs a second order polynomial fit to energy or resolution
// calilbration pairs (gamma or beta)
//
//
void cal_fit_plot(double* eg, double* ec, double* eg_err, int n, double *pec, double *pce, char* opt,int use_err){
int i;
//
//
// Set gamma energy errors to 0
for (i=0; iSetTitle("Energy(keV) vs Channel");
cal_ec->SetTitle("Channel vs Energy(keV)");
cal_ce->GetYaxis()->SetTitle("Channel");
cal_ce->GetXaxis()->SetTitle("Energy (keV)");
cal_ec->GetXaxis()->SetTitle("Channel");
cal_ec->GetYaxis()->SetTitle("Energy (keV)");
cal_ec->SetMarkerStyle(20);
cal_ec->SetMarkerColor(3);
cal_ec->SetMarkerSize(1);
cal_ce->SetMarkerStyle(20);
cal_ce->SetMarkerColor(3);
cal_ce->SetMarkerSize(1);
//
//Perform fits
//
if ( strncmp(opt,"gammaen",7) == 0){
cal_ce->SetTitle("Energy(keV) vs Channel");
cal_ec->SetTitle("Channel vs Energy(keV)");
cal_ce->GetYaxis()->SetTitle("Channel");
cal_ce->GetXaxis()->SetTitle("Energy (keV)");
cal_ec->GetXaxis()->SetTitle("Channel");
cal_ec->GetYaxis()->SetTitle("Energy (keV)");
TCanvas *gecal = new TCanvas("gecal","Gamma energy calibration",500,600);
gecal->Divide(1,2);
gecal->cd(1);
// Draw graphs
cal_ec->Draw("AP");
//Perform pol fit
gStyle->SetOptFit(101);
cal_ec->Fit("pol2");
gPad->Update();
statsbox=(TPaveStats *)(gPad->GetPrimitive("stats"));
statsbox->SetX1NDC(0.15);
statsbox->SetX2NDC(0.4);
statsbox->SetY1NDC(0.5);
statsbox->SetY2NDC(0.85);
statsbox->SetFillColor(41);
statsbox->SetTextFont(2);
statsbox->SetBorderSize(0);
gPad->Draw();
// get fit parameters
TF1 *fitcal = cal_ec->GetFunction("pol2");
pec[0] = fitcal->GetParameter(0);
pec[1] = fitcal->GetParameter(1);
pec[2] = fitcal->GetParameter(2);
pec[3] = fitcal->GetChisquare();
gecal->cd(2);
cal_ce->Draw("AP");
//Perform pol fit
gStyle->SetOptFit(101);
cal_ce->Fit("pol2");
gPad->Update();
statsbox=(TPaveStats *)(gPad->GetPrimitive("stats"));
statsbox->SetX1NDC(0.15);
statsbox->SetX2NDC(0.4);
statsbox->SetY1NDC(0.5);
statsbox->SetY2NDC(0.85);
statsbox->SetFillColor(41);
statsbox->SetTextFont(2);
statsbox->SetBorderSize(0);
gPad->Draw();
// get fit parameters
TF1 *fitcal = cal_ce->GetFunction("pol2");
pce[0] = fitcal->GetParameter(0);
pce[1] = fitcal->GetParameter(1);
pce[2] = fitcal->GetParameter(2);
pce[3] = fitcal->GetChisquare();
}
else if ( strncmp(opt,"betaen",6) == 0){
cal_ce->SetTitle("Energy(keV) vs Channel");
cal_ec->SetTitle("Channel vs Energy(keV)");
cal_ce->GetYaxis()->SetTitle("Channel");
cal_ce->GetXaxis()->SetTitle("Energy (keV)");
cal_ec->GetXaxis()->SetTitle("Channel");
cal_ec->GetYaxis()->SetTitle("Energy (keV)");
TCanvas *becal = new TCanvas("becal","Beta energy calibration",500,600);
becal->Divide(1,2);
becal->cd(1);
// Draw graphs
cal_ec->Draw("AP");
//Perform pol fit
gStyle->SetOptFit(101);
cal_ec->Fit("pol2");
gPad->Update();
statsbox=(TPaveStats *)(gPad->GetPrimitive("stats"));
statsbox->SetX1NDC(0.15);
statsbox->SetX2NDC(0.4);
statsbox->SetY1NDC(0.5);
statsbox->SetY2NDC(0.85);
statsbox->SetFillColor(41);
statsbox->SetTextFont(2);
statsbox->SetBorderSize(0);
gPad->Draw();
// get fit parameters
TF1 *fitcal = cal_ec->GetFunction("pol2");
pec[0] = fitcal->GetParameter(0);
pec[1] = fitcal->GetParameter(1);
pec[2] = fitcal->GetParameter(2);
pec[3] = fitcal->GetChisquare();
becal->cd(2);
cal_ce->Draw("AP");
//Perform pol fit
gStyle->SetOptFit(101);
cal_ce->Fit("pol2");
gPad->Update();
statsbox=(TPaveStats *)(gPad->GetPrimitive("stats"));
statsbox->SetX1NDC(0.15);
statsbox->SetX2NDC(0.4);
statsbox->SetY1NDC(0.5);
statsbox->SetY2NDC(0.85);
statsbox->SetFillColor(41);
statsbox->SetTextFont(2);
statsbox->SetBorderSize(0);
gPad->Draw();
// get fit parameters
TF1 *fitcal = cal_ce->GetFunction("pol2");
pce[0] = fitcal->GetParameter(0);
pce[1] = fitcal->GetParameter(1);
pce[2] = fitcal->GetParameter(2);
pce[3] = fitcal->GetChisquare();
}
else if ( strncmp(opt,"gammares",8) == 0){
cal_ce->SetTitle("Energy(keV) vs FWHM (keV)");
cal_ec->SetTitle("FWHM(keV) vs Energy(keV)");
cal_ce->GetYaxis()->SetTitle("FWHM (keV)");
cal_ce->GetXaxis()->SetTitle("Energy (keV)");
cal_ec->GetXaxis()->SetTitle("FWHM (keV)");
cal_ec->GetYaxis()->SetTitle("Energy (keV)");
TCanvas *grcal = new TCanvas("grcal","Gamma resolution calibration",500,600);
grcal->Divide(1,2);
grcal->cd(1);
// Draw graphs
cal_ec->Draw("AP");
//Perform pol fit
gStyle->SetOptFit(101);
cal_ec->Fit("pol2");
gPad->Update();
statsbox=(TPaveStats *)(gPad->GetPrimitive("stats"));
statsbox->SetX1NDC(0.15);
statsbox->SetX2NDC(0.4);
statsbox->SetY1NDC(0.5);
statsbox->SetY2NDC(0.85);
statsbox->SetFillColor(41);
statsbox->SetTextFont(2);
statsbox->SetBorderSize(0);
gPad->Draw();
// get fit parameters
TF1 *fitcal = cal_ec->GetFunction("pol2");
pec[0] = fitcal->GetParameter(0);
pec[1] = fitcal->GetParameter(1);
pec[2] = fitcal->GetParameter(2);
pec[3] = fitcal->GetChisquare();
grcal->cd(2);
cal_ce->Draw("AP");
//Perform pol fit
gStyle->SetOptFit(101);
cal_ce->Fit("pol2");
gPad->Update();
statsbox=(TPaveStats *)(gPad->GetPrimitive("stats"));
statsbox->SetX1NDC(0.15);
statsbox->SetX2NDC(0.4);
statsbox->SetY1NDC(0.5);
statsbox->SetY2NDC(0.85);
statsbox->SetFillColor(41);
statsbox->SetTextFont(2);
statsbox->SetBorderSize(0);
gPad->Draw();
// get fit parameters
TF1 *fitcal = cal_ce->GetFunction("pol2");
pce[0] = fitcal->GetParameter(0);
pce[1] = fitcal->GetParameter(1);
pce[2] = fitcal->GetParameter(2);
pce[3] = fitcal->GetChisquare();
}
else {
cal_ce->SetTitle("Energy(keV) vs FWHM (keV)");
cal_ec->SetTitle("FWHM(keV) vs Energy(keV)");
cal_ce->GetYaxis()->SetTitle("FWHM (keV)");
cal_ce->GetXaxis()->SetTitle("Energy (keV)");
cal_ec->GetXaxis()->SetTitle("FWHM (keV)");
cal_ec->GetYaxis()->SetTitle("Energy (keV)");
TCanvas *brcal = new TCanvas("brcal","Beta resolution calibration",500,600);
brcal->Divide(1,2);
brcal->cd(1);
// Draw graphs
cal_ec->Draw("AP");
//Perform pol fit
gStyle->SetOptFit(101);
cal_ec->Fit("pol2");
gPad->Update();
statsbox=(TPaveStats *)(gPad->GetPrimitive("stats"));
statsbox->SetX1NDC(0.15);
statsbox->SetX2NDC(0.4);
statsbox->SetY1NDC(0.5);
statsbox->SetY2NDC(0.85);
statsbox->SetFillColor(41);
statsbox->SetTextFont(2);
statsbox->SetBorderSize(0);
gPad->Draw();
// get fit parameters
TF1 *fitcal = cal_ec->GetFunction("pol2");
pec[0] = fitcal->GetParameter(0);
pec[1] = fitcal->GetParameter(1);
pec[2] = fitcal->GetParameter(2);
pec[3] = fitcal->GetChisquare();
brcal->cd(2);
cal_ce->Draw("AP");
//Perform pol fit
gStyle->SetOptFit(101);
cal_ce->Fit("pol2");
gPad->Update();
statsbox=(TPaveStats *)(gPad->GetPrimitive("stats"));
statsbox->SetX1NDC(0.15);
statsbox->SetX2NDC(0.4);
statsbox->SetY1NDC(0.5);
statsbox->SetY2NDC(0.85);
statsbox->SetFillColor(41);
statsbox->SetTextFont(2);
statsbox->SetBorderSize(0);
gPad->Draw();
// get fit parameters
TF1 *fitcal = cal_ce->GetFunction("pol2");
pce[0] = fitcal->GetParameter(0);
pce[1] = fitcal->GetParameter(1);
pce[2] = fitcal->GetParameter(2);
pce[3] = fitcal->GetChisquare();
}
}
// A polynom
//
//
void polconv(double* x, double* y, double * p, int nx){
for (int n = 0; n < nx; n++){
y[n] = p[0]+p[1]*x[n]+p[2]*x[n]**2;
}
}
//
// Converts date and time to number of seconds since 2004-01-01, 00:00.
// Input parameters:
// (char date , char time )
// Return value:
// int
//
int date_time_conv(char date[20], char time[20]){
int year, month,day;
int hour,min,sec;
int dday;
int time_hour,time_min,time_sec;
if (strpbrk(date,"/")!='\0'){
sscanf(date,"%d/%d/%d",&year,&month,&day);
//printf("year: %d month: %d day: %d \n",year,month,day);
}
if (strpbrk(date,"-")!='\0'){
sscanf(date,"%d-%d-%d",&year,&month,&day);
//printf("\n year: %d",year);
}
sscanf(time,"%d:%d:%d",&hour,&min,&sec);
//printf("hour: %d min: %d sec: %d \n",hour,min,sec);
dday=0;
for (int i=2005; i < year;i++){
dday = dday + day_of_year(i,12,31);
}
dday = dday + day_of_year(year,month,day);
//time_hour = dday*24 + hour;
//time_min = dday*24*60 + hour*60 + min;
time_sec = dday*24*60*60 + hour*60*60 + min*60 + sec;
return time_sec;
}
int day_of_year(int year,int month,int day){
// set day of year from month and day
int daytab[2][13] = {
{0,31,28,31,30,31,30,31,31,30,31,30,31},
{0,31,29,31,30,31,30,31,31,30,31,30,31}};
int leap,i;
leap = year%4 == 0 && year%100 !=0 || year%400 == 0;
//printf("\n leap: %d month: %d day: %d",leap,month,day);
for (i=1; i < month; i++){
day += daytab[leap][i];
}
return day;
}
// Performs a second order polynomial fit to energy or resolution
// calilbration pairs (gamma or beta)
//
//
void cal_fit(double* eg, double* ec, double* eg_err, int n, double *pec, double *pce, char* opt,int use_err){
int i;
//
//
// Set gamma energy errors to 0
for (i=0; iSetOptFit(101);
cal_ec->Fit("pol2","Q");
cal_ce->Fit("pol2","Q");
TF1 *fitcal = cal_ec->GetFunction("pol2");
pec[0] = fitcal->GetParameter(0);
pec[1] = fitcal->GetParameter(1);
pec[2] = fitcal->GetParameter(2);
pec[3] = fitcal->GetChisquare();
TF1 *fitcal = cal_ce->GetFunction("pol2");
pce[0] = fitcal->GetParameter(0);
pce[1] = fitcal->GetParameter(1);
pce[2] = fitcal->GetParameter(2);
pce[3] = fitcal->GetChisquare();
}
int get_roi_counts(char* file, int* roi_counts, TCutG *roi[30], char *hname, int* bg_chan_ene, int calplot){
double g_ene[50], g_chan[50], g_err[50]; // Gamma energy-channel callibration pairs with error
int n_gammacal; // Number of gamma energy-channel calibration pairs
double b_ene[50], b_chan[50], b_err[50]; // Beta energy-channel callibration pairs with error
int n_betacal; // Number of beta energy-channel calibration pairs
double roi_b_low[30],roi_b_high[30]; // ROI limits (keV) in beta, corresponding ROI number stored in roi_lim
double roi_g_low[30],roi_g_high[30]; // ROI limits (keV) in gamma
int n_roi; // Number of ROI limits
int roi_lim[30]; // ROI number with limits in roi_b_low, roi_b_high, roi_g_low and roi_g_high
double roi_b_low_c[30], roi_b_high_c[30], roi_g_low_c[30], roi_g_high_c[30];
int croi_b_low_c[30], croi_b_high_c[30], croi_g_low_c[30], croi_g_high_c[30];
double pg_ec[10],pg_ce[10],pb_ec[10],pb_ce[10];
TH2F *hist=(TH2F*)gDirectory->Get(hname);
//Create energy calilbration
n_gammacal=g_energy(file, g_ene, g_chan, g_err);
n_betacal=b_energy(file, b_ene, b_chan, b_err);
if (calplot == 0){
cal_fit(g_ene,g_chan,g_err,n_gammacal,pg_ec,pg_ce,"gammaen",0);
cal_fit(b_ene,b_chan,b_err,n_betacal,pb_ec,pb_ce,"betaen",0);
}
if (calplot == 1){
cal_fit_plot(g_ene,g_chan,g_err,n_gammacal,pg_ec,pg_ce,"gammaen",0);
cal_fit_plot(b_ene,b_chan,b_err,n_betacal,pb_ec,pb_ce,"betaen",0);
}
//get ROI limits in energy
n_roi=ROI_limits(file, roi_lim, roi_b_low, roi_b_high, roi_g_low, roi_g_high);
//printf("ROI limits (keV):\n");
//for (int i=0;iSetPoint(0,croi_b_low_c[i],croi_g_low_c[i]);
roi[i]->SetPoint(1,croi_b_high_c[i],croi_g_low_c[i]);
roi[i]->SetPoint(2,croi_b_high_c[i],croi_g_high_c[i]);
roi[i]->SetPoint(3,croi_b_low_c[i],croi_g_high_c[i]);
roi[i]->SetPoint(4,croi_b_low_c[i],croi_g_low_c[i]);
}
//
//Extract counts in each ROI from the histogram hbg
//
double cont, ccroi[30];
for (i=0;iGetBinContent(x+1,y+1);
for(int i=0;i < n_roi;i++){
if (((x >= croi_b_low_c[i]) && (x<= croi_b_high_c[i])) && ((y >= croi_g_low_c[i]) && (y <= croi_g_high_c[i]))) ccroi[i] += cont; //This gives the same result as the next line
//if (roi[i]->IsInside(x,y)) ccroi[i] +=cont;
}
}
}
for(int i=0;i < n_roi;i++) roi_counts[i]=(int)ccroi[i];
return n_roi;
}
int get_roi_counts_shift(char* file, int* roi_counts, TCutG *roi[30], char *hname, int* bg_chan_ene, double* shift, int calplot){
double g_ene[50], g_chan[50], g_err[50]; // Gamma energy-channel callibration pairs with error
int n_gammacal; // Number of gamma energy-channel calibration pairs
double b_ene[50], b_chan[50], b_err[50]; // Beta energy-channel callibration pairs with error
int n_betacal; // Number of beta energy-channel calibration pairs
double roi_b_low[30],roi_b_high[30]; // ROI limits (keV) in beta, corresponding ROI number stored in roi_lim
double roi_g_low[30],roi_g_high[30]; // ROI limits (keV) in gamma
int n_roi; // Number of ROI limits
int roi_lim[30]; // ROI number with limits in roi_b_low, roi_b_high, roi_g_low and roi_g_high
double roi_b_low_c[30], roi_b_high_c[30], roi_g_low_c[30], roi_g_high_c[30];
int croi_b_low_c[30], croi_b_high_c[30], croi_g_low_c[30], croi_g_high_c[30];
double pg_ec[10],pg_ce[10],pb_ec[10],pb_ce[10];
TH2F *hist=(TH2F*)gDirectory->Get(hname);
//TCanvas *Myc1 = new TCanvas("Myc1",file,800,800);
//hist->Draw();
//Create energy calibration
n_gammacal=g_energy(file, g_ene, g_chan, g_err);
n_betacal=b_energy(file, b_ene, b_chan, b_err);
if (calplot == 0){
cal_fit(g_ene,g_chan,g_err,n_gammacal,pg_ec,pg_ce,"gammaen",0);
cal_fit(b_ene,b_chan,b_err,n_betacal,pb_ec,pb_ce,"betaen",0);
}
if (calplot == 1){
cal_fit_plot(g_ene,g_chan,g_err,n_gammacal,pg_ec,pg_ce,"gammaen",0);
cal_fit_plot(b_ene,b_chan,b_err,n_betacal,pb_ec,pb_ce,"betaen",0);
}
//get ROI limits in energy
n_roi=ROI_limits(file, roi_lim, roi_b_low, roi_b_high, roi_g_low, roi_g_high);
//printf("ROI limits (keV):\n");
//for (int i=0;iSetPoint(0,croi_b_low_c[i],croi_g_low_c[i]);
roi[i]->SetPoint(1,croi_b_high_c[i],croi_g_low_c[i]);
roi[i]->SetPoint(2,croi_b_high_c[i],croi_g_high_c[i]);
roi[i]->SetPoint(3,croi_b_low_c[i],croi_g_high_c[i]);
roi[i]->SetPoint(4,croi_b_low_c[i],croi_g_low_c[i]);
}
//
//Extract counts in each ROI from the histogram hbg
//
double cont, ccroi[30];
for (i=0;iGetBinContent(x+1,y+1);
for(int i=0;i < n_roi;i++){
if (((x >= croi_b_low_c[i]) && (x<= croi_b_high_c[i])) && ((y >= croi_g_low_c[i]) && (y <= croi_g_high_c[i]))) ccroi[i] += cont; //This gives the same result as the next line
//if (roi[i]->IsInside(x,y)) ccroi[i] +=cont;
}
}
}
for(int i=0;i < n_roi;i++) roi_counts[i]=(int)ccroi[i];
return n_roi;
}
int get_own_roi_counts(char* file, int* roi_counts, TCutG *roi[30], char *hname, int* bg_chan_ene, int calplot){
double g_ene[50], g_chan[50], g_err[50]; // Gamma energy-channel callibration pairs with error
int n_gammacal; // Number of gamma energy-channel calibration pairs
double b_ene[50], b_chan[50], b_err[50]; // Beta energy-channel callibration pairs with error
int n_betacal; // Number of beta energy-channel calibration pairs
double pg_ec[10],pg_ce[10],pb_ec[10],pb_ce[10];
int n_roi; // Number of ROI limits
int roi_lim[30]; // ROI number with limits in roi_b_low, roi_b_high, roi_g_low and roi_g_high
double roi_b_low_c[30], roi_b_high_c[30], roi_g_low_c[30], roi_g_high_c[30];
int croi_b_low_c[30], croi_b_high_c[30], croi_g_low_c[30], croi_g_high_c[30];
TH2F *hist=(TH2F*)gDirectory->Get(hname);
n_roi=ROI_limits("param.txt", roi_lim, roi_b_low_c, roi_b_high_c, roi_g_low_c, roi_g_high_c);
//printf("ROI limits (channel user defined):\n");
n_gammacal=g_energy(file, g_ene, g_chan, g_err);
n_betacal=b_energy(file, b_ene, b_chan, b_err);
if (calplot == 1){
cal_fit_plot(g_ene,g_chan,g_err,n_gammacal,pg_ec,pg_ce,"gammaen",0);
cal_fit_plot(b_ene,b_chan,b_err,n_betacal,pb_ec,pb_ce,"betaen",0);
}
for (int i=0;iSetPoint(0,croi_b_low_c[i],croi_g_low_c[i]);
roi[i]->SetPoint(1,croi_b_high_c[i],croi_g_low_c[i]);
roi[i]->SetPoint(2,croi_b_high_c[i],croi_g_high_c[i]);
roi[i]->SetPoint(3,croi_b_low_c[i],croi_g_high_c[i]);
roi[i]->SetPoint(4,croi_b_low_c[i],croi_g_low_c[i]);
}
//
//Extract counts in each ROI from the histogram hbg
//
int cont, ccroi[30];
for (i=0;iGetBinContent(x+1,y+1);
for(int i=0;i < n_roi;i++){
if (((x >= croi_b_low_c[i]) && (x<= croi_b_high_c[i])) && ((y >= croi_g_low_c[i]) && (y <= croi_g_high_c[i]))) ccroi[i] += cont;
}
}
}
for(int i=0;i < n_roi;i++) roi_counts[i]=(int)ccroi[i];
return n_roi;
}
//
//Calculation of gas background subtraction factor
//
double calc_F(double lambda, int tao, double s_live, double s_real, double g_live, double g_real)
{
double f;
f=(g_real/g_live)*(s_live/s_real)*exp(-lambda*tao)*(1-exp(-lambda*s_real))/(1-exp(-lambda*g_real));
return f;
}
double get_xe_counts(int *roi, int *roib, double *R, double sf, double *k, double *Vk, double *LK, int *xrf)
{
//Net counts are stored as follows:
// k[1..10] net counts in ROI 1..10
// Vk[1..10] variance for ROI 1..10
// LK[1..10] critical level for ROI 1..10
//
// Start by calculating net counts in ROI1 (radon):
//
k[1] = roi[0] - sf*roib[0];
Vk[1]= roi[0] + sf*sf*roib[0];
LK[1] = ks*sqrt(Vk[1]-k[1]);
//
//If radon present calculate ROI2 (135Xe) and ROI3 (133 80 keV) correcting for radon
//
if (k[1] > LK[1]){
xrf[0]=1;
k[2] = roi[1] - sf*roib[1]-R[0]*k[1];
Vk[2] = roi[1] + sf*sf*roib[1] + R[0]*R[0]*Vk[1];
k[3] = roi[2] - sf*roib[2] - R[1]*k[1];
Vk[3] = roi[2] + sf*sf*roib[2] + R[1]*R[1]*Vk[1];
LK[3] = ks*sqrt(Vk[3]-k[3]);
// If xe133 present correct also for that in the rest of the ROIs
if (k[3]>LK[3]){
xrf[1]=1;
k[4] = roi[3] - sf*roib[3]-R[2]*k[1]-R[9]*k[3];
Vk[4] = roi[3] + sf*sf*roib[3]+R[2]*R[2]*Vk[1]+R[9]*R[9]*Vk[3];
k[5] = roi[4] - sf*roib[4]-R[3]*k[1]-R[10]*k[3];
Vk[5] = roi[4] + sf*sf*roib[4]+R[3]*R[3]*Vk[1]+R[10]*R[10]*Vk[3];
k[6] = roi[5] - sf*roib[5]-R[4]*k[1]-R[11]*k[3];
Vk[6] = roi[5] + sf*sf*roib[5]+R[4]*R[4]*Vk[1]+R[11]*R[11]*Vk[3];
k[7] = roi[6] - sf*roib[6]-R[5]*k[1]-R[12]*k[3];
Vk[7] = roi[6] + sf*sf*roib[6]+R[5]*R[5]*Vk[1]+R[12]*R[12]*Vk[3];
k[8] = roi[7] - sf*roib[7]-R[6]*k[1]-R[13]*k[3];
Vk[8] = roi[7] + sf*sf*roib[7]+R[6]*R[6]*Vk[1]+R[13]*R[13]*Vk[3];
k[9] = roi[8] - sf*roib[8]-R[7]*k[1]-R[14]*k[3];
Vk[9] = roi[8] + sf*sf*roib[8]+R[7]*R[7]*Vk[1]+R[14]*R[14]*Vk[3];
k[10] = roi[9] - sf*roib[9]-R[8]*k[1]-R[15]*k[3];
Vk[10] = roi[9] + sf*sf*roib[9]+R[8]*R[8]*Vk[1]+R[15]*R[15]*Vk[3];
}
else{
xrf[1]=0;
k[4] = roi[3] - sf*roib[3]-R[2]*k[1];
Vk[4] = roi[3] + sf*sf*roib[3]+R[2]*R[2]*Vk[1];
k[5] = roi[4] - sf*roib[4]-R[3]*k[1];
Vk[5] = roi[4] + sf*sf*roib[4]+R[3]*R[3]*Vk[1];
k[6] = roi[5] - sf*roib[5]-R[4]*k[1];
Vk[6] = roi[5] + sf*sf*roib[5]+R[4]*R[4]*Vk[1];
k[7] = roi[6] - sf*roib[6]-R[5]*k[1];
Vk[7] = roi[6] + sf*sf*roib[6]+R[5]*R[5]*Vk[1];
k[8] = roi[7] - sf*roib[7]-R[6]*k[1];
Vk[8] = roi[7] + sf*sf*roib[7]+R[6]*R[6]*Vk[1];
k[9] = roi[8] - sf*roib[8]-R[7]*k[1];
Vk[9] = roi[8] + sf*sf*roib[8]+R[7]*R[7]*Vk[1];
k[10] = roi[9] - sf*roib[9]-R[8]*k[1];
Vk[10] = roi[9] + sf*sf*roib[9]+R[8]*R[8]*Vk[1];
}
}
// If no radon present:
else{
k[2] = roi[1] - sf*roib[1];
Vk[2] = roi[1] + sf*sf*roib[1];
k[3] = roi[2] - sf*roib[2];
Vk[3] = roi[2] + sf*sf*roib[2];
LK[3] = ks*sqrt(Vk[3]-k[3]);
xrf[0]=0;
// If xe133 present correct also for that in the rest of the ROIs
if (k[3]>LK[3]){
xrf[1]=1;
k[4] = roi[3] - sf*roib[3]-R[9]*k[3];
Vk[4] = roi[3] + sf*sf*roib[3]+R[9]*R[9]*Vk[3];
k[5] = roi[4] - sf*roib[4]-R[10]*k[3];
Vk[5] = roi[4] + sf*sf*roib[4]+R[10]*R[10]*Vk[3];
k[6] = roi[5] - sf*roib[5]-R[11]*k[3];
Vk[6] = roi[5] + sf*sf*roib[5]+R[11]*R[11]*Vk[3];
k[7] = roi[6] - sf*roib[6]-R[12]*k[3];
Vk[7] = roi[6] + sf*sf*roib[6]+R[12]*R[12]*Vk[3];
k[8] = roi[7] - sf*roib[7]-R[13]*k[3];
Vk[8] = roi[7] + sf*sf*roib[7]+R[13]*R[13]*Vk[3];
k[9] = roi[8] - sf*roib[8]-R[14]*k[3];
Vk[9] = roi[8] + sf*sf*roib[8]+R[14]*R[14]*Vk[3];
k[10] = roi[9] - sf*roib[9]-R[15]*k[3];
Vk[10] = roi[9] + sf*sf*roib[9]+R[15]*R[15]*Vk[3];
}
else{
xrf[1]=0;
k[4] = roi[3] - sf*roib[3];
Vk[4] = roi[3] + sf*sf*roib[3];
k[5] = roi[4] - sf*roib[4];
Vk[5] = roi[4] + sf*sf*roib[4];
k[6] = roi[5] - sf*roib[5];
Vk[6] = roi[5] + sf*sf*roib[5];
k[7] = roi[6] - sf*roib[6];
Vk[7] = roi[6] + sf*sf*roib[6];
k[8] = roi[7] - sf*roib[7];
Vk[8] = roi[7] + sf*sf*roib[7];
k[9] = roi[8] - sf*roib[8];
Vk[9] = roi[8] + sf*sf*roib[8];
k[10] = roi[9] - sf*roib[9];
Vk[10] = roi[9] + sf*sf*roib[9];
}
}
LK[2]= ks*sqrt(Vk[2]-k[2]);
LK[3]= ks*sqrt(Vk[3]-k[3]);
LK[4]= ks*sqrt(Vk[4]-k[4]);
LK[5]= ks*sqrt(Vk[5]-k[5]);
LK[6]= ks*sqrt(Vk[6]-k[6]);
LK[7]= ks*sqrt(Vk[7]-k[7]);
LK[8]= ks*sqrt(Vk[8]-k[8]);
LK[9]= ks*sqrt(Vk[9]-k[9]);
LK[10]= ks*sqrt(Vk[10]-k[10]);
}
double get_net_counts(double *ks, double *kg, double *Vks, double *Vkg, double *nc, double *Vn, double *LCn, double *F)
{
// ROI 1
nc[1] = ks[1] - F[4]*kg[1];
Vn[1] = Vks[1]+ F[4]*F[4]*Vkg[1];
LCn[1] = Vn[1] - nc[1];
// ROI 2
nc[2] = ks[2] - F[1]*kg[2];
Vn[2] = Vks[2]+ F[1]*F[1]*Vkg[2];
LCn[2] = Vn[2] - nc[2];
// ROI 3
nc[3] = ks[3] - F[0]*kg[3];
Vn[3] = Vks[3]+ F[0]*F[0]*Vkg[3];
LCn[3] = Vn[3] - nc[3];
// ROI 4
nc[4] = ks[4] - F[0]*kg[4];
Vn[4] = Vks[4]+ F[0]*F[0]*Vkg[4];
LCn[4] = Vn[4] - nc[4];
// ROI 5
nc[5] = ks[5] - F[2]*kg[5];
Vn[5] = Vks[5]+ F[2]*F[2]*Vkg[5];
LCn[5] = Vn[5] - nc[5];
// ROI 6
nc[6] = ks[6] - F[3]*kg[6];
Vn[6] = Vks[6]+ F[3]*F[3]*Vkg[6];
LCn[6] = Vn[6] - nc[6];
// ROI 7
nc[7] = ks[7] - F[0]*kg[7];
Vn[7] = Vks[7]+ F[0]*F[0]*Vkg[7];
LCn[7] = Vn[7] - nc[7];
// ROI 8
nc[8] = ks[8] - F[0]*kg[8];
Vn[8] = Vks[8]+ F[0]*F[0]*Vkg[8];
LCn[8] = Vn[8] - nc[8];
// ROI 9
nc[9] = ks[9] - F[0]*kg[9];
Vn[9] = Vks[9]+ F[0]*F[0]*Vkg[9];
LCn[9] = Vn[9] - nc[9];
// ROI 10
nc[10] = ks[10] - F[0]*kg[10];
Vn[10] = Vks[10]+ F[0]*F[0]*Vkg[10];
LCn[10] = Vn[10] - nc[10];
}
int get_conc_roi(double *nc, double *Vn, double *eff, int tc, int tp, int ta, double *Vxe, double *croi, double *cerr_roi, double *kroi)
{
croi[1] = (nc[1]/(eff[0]*bg_br_222Rn))*calc_conc(Rn222_lam,tc,tp,ta,Vxe);
//printf("nc1: %d\n",nc[1]);
//printf("nc3: %d\n",nc[3]);
cerr_roi[1]=0;
//if (nc[1]!=0){
//cerr_roi[1] = croi[1]*sqrt(Vn[1])/nc[1];
//}
croi[2] = (nc[2]/(eff[0]*bg_br_135Xe))*calc_conc(Xe135_lam,tc,tp,ta,Vxe);
cerr_roi[2] = croi[2]*sqrt(Vn[2])/nc[2];
croi[3] = (nc[3]/(eff[1]*bg_br_133Xe_80))*calc_conc(Xe133_lam,tc,tp,ta,Vxe);
cerr_roi[3] = croi[3]*sqrt(Vn[3])/nc[3];
croi[4] = (nc[4]/(eff[2]*bg_br_133Xe_30))*calc_conc(Xe133_lam,tc,tp,ta,Vxe);
cerr_roi[4] = croi[4]*sqrt(Vn[4])/nc[4];
croi[5] = (nc[5]/(eff[3]*bg_br_131mXe_30))*calc_conc(Xe131m_lam,tc,tp,ta,Vxe);
cerr_roi[5] = croi[5]*sqrt(Vn[5])/nc[5];
croi[6] = (nc[6]/(eff[4]*bg_br_133mXe_30))*calc_conc(Xe133m_lam,tc,tp,ta,Vxe);
cerr_roi[6] = croi[6]*sqrt(Vn[6])/nc[6];
croi[7] = (nc[7]/(eff[5]*bg_br_133Xe_30))*calc_conc(Xe133_lam,tc,tp,ta,Vxe);
cerr_roi[7] = croi[7]*sqrt(Vn[7])/nc[7];
croi[8] = (nc[8]/(eff[6]*bg_br_133Xe_30))*calc_conc(Xe133_lam,tc,tp,ta,Vxe);
cerr_roi[8] = croi[8]*sqrt(Vn[8])/nc[8];
croi[9] = (nc[9]/(eff[7]*bg_br_133Xe_30))*calc_conc(Xe133_lam,tc,tp,ta,Vxe);
cerr_roi[9] = croi[9]*sqrt(Vn[9])/nc[9];
croi[10] = (nc[10]/(eff[8]*bg_br_133Xe_30))*calc_conc(Xe133_lam,tc,tp,ta,Vxe);
cerr_roi[10] = croi[10]*sqrt(Vn[10])/nc[10];
kroi[1]=0;
//kroi[1] = nc[1]/croi[1];
kroi[2] = nc[2]/croi[2];
kroi[3] = nc[3]/croi[3];
kroi[4] = nc[4]/croi[4];
kroi[5] = nc[5]/croi[5];
kroi[6] = nc[6]/croi[6];
kroi[7] = nc[7]/croi[7];
kroi[8] = nc[8]/croi[8];
kroi[9] = nc[9]/croi[9];
kroi[10] = nc[10]/croi[10];
//Investigate if metastables present
//calulate critical limits
double lc3 = ks*sqrt(Vn[3]- nc[3]);
double lc4 = ks*sqrt(Vn[4]-nc[4]);
double lc5 = ks*sqrt(Vn[5]-nc[5]);
double lc6 = ks*sqrt(Vn[6]-nc[6]);
double lc7 = ks*sqrt(Vn[7]-nc[7]);
double lc8 = ks*sqrt(Vn[8]-nc[8]);
double lc9 = ks*sqrt(Vn[9]-nc[9]);
double lc10 = ks*sqrt(Vn[10]-nc[10]);
int cas;
if (lc5 < nc[5]){
if (lc6 < nc[6]){
//xe131m and xe133m detected
cas = 1;
}
else {
//Xe131m only detected
cas = 2;
}
}
else if (lc6 < nc[6]){
//Xe133m only detected
cas = 3;
}
else {
//no metastables detected
cas = 4;
}
//printf("cas: %d\n",cas);
double xif = xi(Xe133m_lam,Xe133_lam,tc,tp,ta);
// printf("xi: %f \n", xif);
croi[3] = croi[3] - xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*croi[6];
cerr_roi[3] = sqrt(cerr_roi[3]*cerr_roi[3] + xif*xif*cerr_roi[6]*cerr_roi[6]);
croi[4] = croi[4] - xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*croi[6];
cerr_roi[4] = sqrt(cerr_roi[4]*cerr_roi[4] + xif*xif*cerr_roi[6]*cerr_roi[6]);
croi[7] = croi[7] - xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*croi[6];
cerr_roi[7] = sqrt(cerr_roi[7]*cerr_roi[7] + xif*xif*cerr_roi[6]*cerr_roi[6]);
croi[8] = croi[8] - xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*croi[6];
cerr_roi[8] = sqrt(cerr_roi[8]*cerr_roi[8] + xif*xif*cerr_roi[6]*cerr_roi[6]);
croi[9] = croi[9] - xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*croi[6];
cerr_roi[9] = sqrt(cerr_roi[9]*cerr_roi[9] + xif*xif*cerr_roi[6]*cerr_roi[6]);
croi[10] = croi[10] - xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*croi[6];
cerr_roi[10] = sqrt(cerr_roi[10]*cerr_roi[10] + xif*xif*cerr_roi[6]*cerr_roi[6]);
return cas;
}
double calc_conc(double lambda, int tc, int tp, int ta, double *Vxe){
double Xe_air=0.087; //Volume (ml) of Xe per m3 of air
double cf;
cf = 1000*lambda*lambda/((1-exp(-lambda*tc))*exp(-lambda*tp)*(1-exp(-lambda*ta)))*(tc/(Vxe[0]/Xe_air));
return cf;
}
double xi(double lambda1, double lambda2, int tc, int tp, int ta){
double fc1=1-exp(-lambda1*tc); double fp1 = exp(-lambda1*tp); double fa1 = 1-exp(-lambda1*ta);
double fc2=1-exp(-lambda2*tc); double fp2 = exp(-lambda2*tp); double fa2 = 1-exp(-lambda2*ta);
return lambda2*lambda2/lambda1/lambda1/(lambda2-lambda1)*fc1/fc2*(lambda2*fp1*fa1/fp2/fa2 - lambda1);
}
void get_conc(int cas, double *croi, double *cerr_roi, double *conc, int tc, int tp, int ta)
//void get_conc(int cas, double *croi, double *cerr_roi, double *k_s, double *LK_s, double *k_g, double *LK_g, double *conc, int tc, int tp, int ta)
{
//
//Returns isotopic xenon concentrations with errors:
//
// conc[0]: 133 conc
// conc[1]: 133 error
// conc[2]: 131m conc
// conc[3]: 131m error
// conc[4]: 133m conc
// conc[5]: 133m error
// conc[6]: 135 conc
// conc[7]: 135 error
//
//
//
//
double C_135 = croi[2];
double dC_135 = cerr_roi[2];
double C_131m = croi[5];
double dC_131m = cerr_roi[5];
double C_133m = croi[6];
double dC_133m = cerr_roi[6];
if (cas ==1){
//Xe131m and xe133m present
double varsum = 1/(cerr_roi[3]*cerr_roi[3])+ 1/(cerr_roi[7]*cerr_roi[7])+1/(cerr_roi[8]*cerr_roi[8]);
double C_133 = (croi[3]/(cerr_roi[3]*cerr_roi[3])+ croi[7]/(cerr_roi[7]*cerr_roi[7]) + croi[8]/(cerr_roi[8]*cerr_roi[8]))/varsum;
}
if (cas == 2){
//Xe131m only present
double varsum = 1/(cerr_roi[3]*cerr_roi[3])+ 1/(cerr_roi[7]*cerr_roi[7])+1/(cerr_roi[9]*cerr_roi[9]);
double C_133 = (croi[3]/(cerr_roi[3]*cerr_roi[3])+ croi[7]/(cerr_roi[7]*cerr_roi[7]) + croi[9]/(cerr_roi[9]*cerr_roi[9]))/varsum;
}
if (cas == 3){
//Xe133m only present
double varsum = 1/(cerr_roi[3]*cerr_roi[3])+ 1/(cerr_roi[8]*cerr_roi[8])+1/(cerr_roi[10]*cerr_roi[10]);
double C_133 = (croi[3]/(cerr_roi[3]*cerr_roi[3])+ croi[8]/(cerr_roi[8]*cerr_roi[8]) + croi[10]/(cerr_roi[10]*cerr_roi[10]))/varsum;
C_133 = C_133 - xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*C_133m;
}
if (cas ==4){
//No metastables present
double varsum = 1/(cerr_roi[3]*cerr_roi[3])+ 1/(cerr_roi[4]*cerr_roi[4]);
double C_133 = (croi[3]/(cerr_roi[3]*cerr_roi[3])+ croi[4]/(cerr_roi[4]*cerr_roi[4]))/varsum;
}
double dC_133 = sqrt(1/varsum);
conc[0]=C_133; conc[1]= dC_133; conc[2]=C_131m; conc[3]=dC_131m;
conc[4]=C_133m; conc[5]=dC_133m; conc[6]=C_135; conc[7]=dC_135;
}
void get_LC_n(double *nc, double *Vn, double *kroi, int cas, double *lc, int tc, int tp, int ta){
double ks = 1.6449;
//
//Returns LC:
// lc[0] : 133Xe
// lc[1] : 131mXe
// lc[2] : 133mXe
// lc[3] : 135Xe
//
double LC_135 = ks*sqrt(Vn[2]-nc[2])/kroi[2];
double LC_131m = ks*sqrt(Vn[5]-nc[5])/kroi[5];
double LC_133m = ks*sqrt(Vn[6]-nc[6])/kroi[6];
double a3 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[3]/kroi[6];
double a4 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[4]/kroi[6];
double a7 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[7]/kroi[6];
double a8 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[8]/kroi[6];
double a9 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[9]/kroi[6];
double a10 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[10]/kroi[6];
double lc3 = ks*sqrt(Vn[3]+a3*a3*Vn[6]-nc[3]+a3*nc[6])/kroi[3];
double lc4 = ks*sqrt(Vn[4]+a4*a4*Vn[6]-nc[4]+a4*nc[6])/kroi[4];
double lc7 = ks*sqrt(Vn[7]+a7*a7*Vn[6]-nc[7]+a7*nc[6])/kroi[7];
double lc8 = ks*sqrt(Vn[8]+a8*a8*Vn[6]-nc[8]+a8*nc[6])/kroi[8];
double lc9 = ks*sqrt(Vn[9]+a9*a9*Vn[6]-nc[9]+a9*nc[6])/kroi[9];
double lc10 = ks*sqrt(Vn[10]+a10*a10*Vn[6]-nc[10]+a10*nc[6])/kroi[10];
if (cas == 4){
//Use ROI 3 and 4, no metastables present
double LC_133 = 1/sqrt(1/lc3/lc3 + 1/lc4/lc4);
}
// Use ROI 3,7 and 9, 131m present
else if (cas == 2){
double LC_133 = 1/sqrt(1/lc3/lc3 + 1/lc7/lc7 + 1/lc9/lc9);
}
// Use ROI 3,8 and 10, 133m present
else if (cas == 3){
double LC_133 = 1/sqrt(1/lc3/lc3 + 1/lc8/lc8 + 1/lc10/lc10);
}
// Use ROI 3,7 and 8, 133m and 131m present
else if (cas == 1){
LC_133 = 1/sqrt( 1/lc3/lc3 + 1/lc7/lc7 + 1/lc8/lc8);
}
lc[0]=LC_133; lc[1]=LC_131m; lc[2]=LC_133m; lc[3]=LC_135;
}
void get_LCA_n(double *nc, double *Vn, double *karoi, int cas, double *lca){
double ks = 1.6449;
//
//Returns LCA:
// lc[0] : 133Xe
// lc[1] : 131mXe
// lc[2] : 133mXe
// lc[3] : 135Xe
// lc[4] : 222Rn in equilibrium
//
double LCA_222 = ks*sqrt(Vn[1]-nc[1])/karoi[1];
double LCA_135 = ks*sqrt(Vn[2]-nc[2])/karoi[2];
double LCA_131m = ks*sqrt(Vn[5]-nc[5])/karoi[5];
double LCA_133m = ks*sqrt(Vn[6]-nc[6])/karoi[6];
// double a3 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[3]/kroi[6];
//double a4 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[4]/kroi[6];
//double a7 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[7]/kroi[6];
//double a8 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[8]/kroi[6];
//double a9 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[9]/kroi[6];
//double a10 = xi(Xe133m_lam,Xe133_lam,tc,tp,ta)*kroi[10]/kroi[6];
//double lc3 = ks*sqrt(Vn[3]+a3*a3*Vn[6]-nc[3]+a3*nc[6])/kroi[3];
//double lc4 = ks*sqrt(Vn[4]+a4*a4*Vn[6]-nc[4]+a4*nc[6])/kroi[4];
//double lc7 = ks*sqrt(Vn[7]+a7*a7*Vn[6]-nc[7]+a7*nc[6])/kroi[7];
//double lc8 = ks*sqrt(Vn[8]+a8*a8*Vn[6]-nc[8]+a8*nc[6])/kroi[8];
//double lc9 = ks*sqrt(Vn[9]+a9*a9*Vn[6]-nc[9]+a9*nc[6])/kroi[9];
//double lc10 = ks*sqrt(Vn[10]+a10*a10*Vn[6]-nc[10]+a10*nc[6])/kroi[10];
double lca3 = ks*sqrt(Vn[3]-nc[3])/karoi[3];
double lca4 = ks*sqrt(Vn[4]-nc[4])/karoi[4];
double lca7 = ks*sqrt(Vn[7]-nc[7])/karoi[7];
double lca8 = ks*sqrt(Vn[8]-nc[8])/karoi[8];
double lca9 = ks*sqrt(Vn[9]-nc[9])/karoi[9];
double lca10 = ks*sqrt(Vn[10]-nc[10])/karoi[10];
if (cas == 4){
//Use ROI 3 and 4, no metastables present
double LCA_133 = 1/sqrt(1/lca3/lca3 + 1/lca4/lca4);
}
// Use ROI 3,7 and 9, 131m present
else if (cas == 2){
double LCA_133 = 1/sqrt(1/lca3/lca3 + 1/lca7/lca7 + 1/lca9/lca9);
}
// Use ROI 3,8 and 10, 133m present
else if (cas == 3){
double LCA_133 = 1/sqrt(1/lca3/lca3 + 1/lca8/lca8 + 1/lca10/lca10);
}
// Use ROI 3,7 and 8, 133m and 131m present
else if (cas == 1){
LCA_133 = 1/sqrt( 1/lca3/lca3 + 1/lca7/lca7 + 1/lca8/lca8);
}
lca[0]=LCA_133; lca[1]=LCA_131m; lca[2]=LCA_133m; lca[3]=LCA_135; lca[4]=LCA_222;
}
void get_LC(double *LCn, double *kroi, int cas, double *lc){
double ks = 1.6449;
//
//Returns LC:
// lc[0] : 133Xe
// lc[1] : 131mXe
// lc[2] : 133mXe
// lc[3] : 135Xe
//
double LC_135 = ks*sqrt(Vn[2]-nc[2])/kroi[2];
double LC_131m = ks*sqrt(Vn[5]-nc[5])/kroi[5];
double LC_133m = ks*sqrt(Vn[6]-nc[6])/kroi[6];
if (cas == 4){
//Use ROI 3 and 4, no metastables present
double LC_133 = ks/sqrt((kroi[3]*kroi[3])/LCn[3]+ (kroi[4]*kroi[4])/LCn[4]);
}
// Use ROI 3,7 and 9, 131m present
else if (cas == 2){
double LC_133 = ks/sqrt((kroi[3]*kroi[3])/LCn[3]+ (kroi[7]*kroi[7])/LCn[7]+ (kroi[9]*kroi[9])/LCn[9]);
}
// Use ROI 3,8 and 10, 133m present
else if (cas == 3){
double LC_133 = ks/sqrt((kroi[3]*kroi[3])/LCn[3]+ (kroi[8]*kroi[8])/LCn[8]+ (kroi[10]*kroi[10])/LCn[10]);
}
// Use ROI 3,7 and 8, 133m and 131m present
else if (cas == 1){
double LC_133 = ks/sqrt((kroi[3]*kroi[3])/LCn[3]+ (kroi[7]*kroi[7])/LCn[7]+ (kroi[8]*kroi[8])/LCn[8]);
}
lc[0]=LC_133; lc[1]=LC_131m; lc[2]=LC_133m; lc[3]=LC_135;
}
void get_MDC(double *LCn, double *lc, double *kroi, int cas, double *mdc){
//
//Returns MDC:
// mdc[0] : 133Xe
// mdc[1] : 131mXe
// mdc[2] : 133mXe
// mdc[3] : 135Xe
//
double MDC_133;
double MDC_135 = lc[3]+ks*sqrt((LCn[2]+LD(LCn[2]))/(kroi[2]*kroi[2]));
double MDC_131m = lc[1]+ks*sqrt((LCn[5]+LD(LCn[5]))/(kroi[5]*kroi[5]));
double MDC_133m = lc[2]+ks*sqrt((LCn[6]+LD(LCn[6]))/(kroi[6]*kroi[6]));
if (cas == 4){
//Use ROI 3 and 4, no metastables present
double MDC_133 = lc[0]+ks*sqrt(
((LCn[3]+LD(LCn[3]))*(LCn[4]+LD(LCn[4])))/
((kroi[3]*kroi[3])*(LCn[4]+LD(LCn[4]))+
(kroi[4]*kroi[4])*(LCn[3]+LD(LCn[3])))
);
}
// Use ROI 3,7 and 9, 131m present
else if (cas == 2){
double vd3 = (LCn[3]+LD(LCn[3]))/(kroi[3]*kroi[3]);
double vd7_9 = (LCn[7]+LD(LCn[7]))*(LCn[9]+LD(LCn[9]))/
((kroi[7]*kroi[7])*(LCn[9]+LD(LCn[9]))+
(kroi[9]*kroi[9])*(LCn[7]+LD(LCn[7])));
double vd3_7_9=1/(1/vd3 + 1/vd7_9);
MDC_133 = ks/sqrt((kroi[3]*kroi[3])/LCn[3]+
(kroi[7]*kroi[7])/LCn[7]+
(kroi[9]*kroi[9])/LCn[9])+
ks*sqrt(vd3_7_9);
}
// Use ROI 3,8 and 10, 133m present
else if (cas == 3){
double vd3 = (LCn[3]+LD(LCn[3]))/(kroi[3]*kroi[3]);
double vd8_10 = (LCn[8]+LD(LCn[8]))*(LCn[10]+LD(LCn[10]))/
((kroi[8]*kroi[8])*(LCn[10]+LD(LCn[10]))+
(kroi[10]*kroi[10])*(LCn[8]+LD(LCn[8])));
double vd3_8_10=1/(1/vd3 + 1/vd8_10);
MDC_133 = ks/sqrt((kroi[3]*kroi[3])/LCn[3]+
(kroi[8]*kroi[8])/LCn[8]+
(kroi[10]*kroi[10])/LCn[10])+
ks*sqrt(vd3_8_10);
}
// Use ROI 3,7 and 8, 133m and 131m present
else if (cas == 1){
double vd3 = (LCn[3]+LD(LCn[3]))/(kroi[3]*kroi[3]);
double vd7_8 = (LCn[7]+LD(LCn[7]))*(LCn[8]+LD(LCn[8]))/
((kroi[7]*kroi[7])*(LCn[8]+LD(LCn[8]))+
(kroi[8]*kroi[8])*(LCn[7]+LD(LCn[7])));
double vd3_7_8=1/(1/vd3 + 1/vd7_8);
MDC_133 = ks/sqrt((kroi[3]*kroi[3])/LCn[3]+
(kroi[7]*kroi[7])/LCn[7]+
(kroi[8]*kroi[8])/LCn[8])+
ks*sqrt(vd3_7_8);
}
mdc[0]=MDC_133; mdc[1]=MDC_131m; mdc[2]=MDC_133m; mdc[3]=MDC_135;
}
void get_MDA(double *LCn, double *lca, double *karoi, int cas, double *mda){
//
//Returns MDA:
// mda[0] : 133Xe
// mda[1] : 131mXe
// mda[2] : 133mXe
// mda[3] : 135Xe
// mda[4] : 222Rn in equilibrium
//
double MDA_133;
double MDA_135 = lca[3]+ks*sqrt((LCn[2]+LD(LCn[2]))/(karoi[2]*karoi[2]));
double MDA_131m = lca[1]+ks*sqrt((LCn[5]+LD(LCn[5]))/(karoi[5]*karoi[5]));
double MDA_133m = lca[2]+ks*sqrt((LCn[6]+LD(LCn[6]))/(karoi[6]*karoi[6]));
double MDA_222 = lca[1]+ks*sqrt((LCn[1]+LD(LCn[1]))/(karoi[1]*karoi[1]));
if (cas == 4){
//Use ROI 3 and 4, no metastables present
double MDA_133 = lca[0]+ks*sqrt(
((LCn[3]+LD(LCn[3]))*(LCn[4]+LD(LCn[4])))/
((karoi[3]*karoi[3])*(LCn[4]+LD(LCn[4]))+
(karoi[4]*karoi[4])*(LCn[3]+LD(LCn[3])))
);
}
// Use ROI 3,7 and 9, 131m present
else if (cas == 2){
double vd3 = (LCn[3]+LD(LCn[3]))/(karoi[3]*karoi[3]);
double vd7_9 = (LCn[7]+LD(LCn[7]))*(LCn[9]+LD(LCn[9]))/
((karoi[7]*karoi[7])*(LCn[9]+LD(LCn[9]))+
(karoi[9]*karoi[9])*(LCn[7]+LD(LCn[7])));
double vd3_7_9=1/(1/vd3 + 1/vd7_9);
MDA_133 = ks/sqrt((karoi[3]*karoi[3])/LCn[3]+
(karoi[7]*karoi[7])/LCn[7]+
(karoi[9]*karoi[9])/LCn[9])+
ks*sqrt(vd3_7_9);
}
// Use ROI 3,8 and 10, 133m present
else if (cas == 3){
double vd3 = (LCn[3]+LD(LCn[3]))/(karoi[3]*karoi[3]);
double vd8_10 = (LCn[8]+LD(LCn[8]))*(LCn[10]+LD(LCn[10]))/
((karoi[8]*karoi[8])*(LCn[10]+LD(LCn[10]))+
(karoi[10]*karoi[10])*(LCn[8]+LD(LCn[8])));
double vd3_8_10=1/(1/vd3 + 1/vd8_10);
MDA_133 = ks/sqrt((karoi[3]*karoi[3])/LCn[3]+
(karoi[8]*karoi[8])/LCn[8]+
(karoi[10]*karoi[10])/LCn[10])+
ks*sqrt(vd3_8_10);
}
// Use ROI 3,7 and 8, 133m and 131m present
else if (cas == 1){
double vd3 = (LCn[3]+LD(LCn[3]))/(karoi[3]*karoi[3]);
double vd7_8 = (LCn[7]+LD(LCn[7]))*(LCn[8]+LD(LCn[8]))/
((karoi[7]*karoi[7])*(LCn[8]+LD(LCn[8]))+
(karoi[8]*karoi[8])*(LCn[7]+LD(LCn[7])));
double vd3_7_8=1/(1/vd3 + 1/vd7_8);
MDA_133 = ks/sqrt((karoi[3]*karoi[3])/LCn[3]+
(karoi[7]*karoi[7])/LCn[7]+
(karoi[8]*karoi[8])/LCn[8])+
ks*sqrt(vd3_7_8);
}
mda[0]=MDA_133; mda[1]=MDA_131m; mda[2]=MDA_133m; mda[3]=MDA_135; mda[4]=MDA_222;
}
double LD(double x){
return ks*ks+2*ks*sqrt(abs(x));
}
void get_act_roi(double *nc, double *Vn, double *eff, int ta, double *aroi, double *aerr_roi, double *karoi)
{
aroi[1] = (nc[1]/(eff[0]*bg_br_222Rn))*calc_act(Rn222_lam,ta);
aerr_roi[1] = aroi[1]*sqrt(Vn[1])/nc[1];
aroi[2] = (nc[2]/(eff[0]*bg_br_135Xe))*calc_act(Xe135_lam,ta);
aerr_roi[2] = aroi[2]*sqrt(Vn[2])/nc[2];
aroi[3] = (nc[3]/(eff[1]*bg_br_133Xe_80))*calc_act(Xe133_lam,ta);
aerr_roi[3] = aroi[3]*sqrt(Vn[3])/nc[3];
aroi[4] = (nc[4]/(eff[2]*bg_br_133Xe_30))*calc_act(Xe133_lam,ta);
aerr_roi[4] = aroi[4]*sqrt(Vn[4])/nc[4];
aroi[5] = (nc[5]/(eff[3]*bg_br_131mXe_30))*calc_act(Xe131m_lam,ta);
aerr_roi[5] = aroi[5]*sqrt(Vn[5])/nc[5];
aroi[6] = (nc[6]/(eff[4]*bg_br_133mXe_30))*calc_act(Xe133m_lam,ta);
aerr_roi[6] = aroi[6]*sqrt(Vn[6])/nc[6];
aroi[7] = (nc[7]/(eff[5]*bg_br_133Xe_30))*calc_act(Xe133_lam,ta);
aerr_roi[7] = aroi[7]*sqrt(Vn[7])/nc[7];
aroi[8] = (nc[8]/(eff[6]*bg_br_133Xe_30))*calc_act(Xe133_lam,ta);
aerr_roi[8] = aroi[8]*sqrt(Vn[8])/nc[8];
aroi[9] = (nc[9]/(eff[7]*bg_br_133Xe_30))*calc_act(Xe133_lam,ta);
aerr_roi[9] = aroi[9]*sqrt(Vn[9])/nc[9];
aroi[10] = (nc[10]/(eff[8]*bg_br_133Xe_30))*calc_act(Xe133_lam,ta);
aerr_roi[10] = aroi[10]*sqrt(Vn[10])/nc[10];
karoi[1] = nc[1]/aroi[1];
karoi[2] = nc[2]/aroi[2];
karoi[3] = nc[3]/aroi[3];
karoi[4] = nc[4]/aroi[4];
karoi[5] = nc[5]/aroi[5];
karoi[6] = nc[6]/aroi[6];
karoi[7] = nc[7]/aroi[7];
karoi[8] = nc[8]/aroi[8];
karoi[9] = nc[9]/aroi[9];
karoi[10] = nc[10]/aroi[10];
}
double calc_act(double lambda, int ta){
double af = 1000*lambda/(1-exp(-lambda*ta));
return af;
}
double fitflin(double *, double *);
double fitflinstep(double *, double *);
double fitfquad(double *, double *);
double peak(double *, double *);
double background(double *, double *);
double backgroundstep(double *, double *);
double background2(double *, double *);
void fitpeak(char histname[80], double x1, double x2,
double scale, double mean, double sigma,int lin){
double slope,bgr,bgrstart,bgrslope,centroid,fwhm,
dummy,quadratic,step,area;
int i,j;
TH1F *hist=(TH1F*)gDirectory->Get(histname);
TCanvas *FitPlot = new TCanvas("FitPlot","fitpeak",1);
FitPlot->SetFillColor(0);
hist->Draw();
if(lin==2){
TF1 *f1 = new TF1("f1",fitfquad,x1,x2,8);
}
else{
if(lin==1){
TF1 *f1 = new TF1("f1",fitflin,x1,x2,8);
}
else if(lin==3){
TF1 *f1 = new TF1("f1",fitflinstep,x1,x2,8);
}
else{
printf("\nError in last parameter: should be 1, 2 or 3. Exiting...\n");
break;
}
}
f1->SetParNames("bgrstart","bgrslope","area","centroid",
"FWHM","","quadratic","step ht");
slope=(hist->GetBinContent(x2)-hist->GetBinContent(x1))/(x2-x1);
printf("\nBackground: at xstart %d, at xend %d, slope %.3f\n\n",
hist->GetBinContent(x1),
hist->GetBinContent(x2),
slope);
f1->SetParameter(0,hist->GetBinContent(x1));
f1->SetParameter(1,slope);
f1->SetParameter(2,scale);
f1->SetParameter(3,mean);
f1->SetParameter(4,sigma);
f1->SetParameter(5,x1);
f1->SetParameter(6,0);
f1->SetParameter(7,0);
if(lin==1||lin==2){
f1->FixParameter(7,0);//step fixed to zero unless linear+step bgr specified
}
if(lin==1||lin==3){
f1->FixParameter(6,0);//quadratic term fixed to zero for linear and linear+step bgr
}
f1->FixParameter(5,x1);
f1->SetParLimits(0,0,100*hist->GetBinContent(x1));
hist->Fit("f1","R");
bgrstart=f1->GetParameter(0);
bgrslope=f1->GetParameter(1);
area=f1->GetParameter(2);
centroid=f1->GetParameter(3);
fwhm=f1->GetParameter(4);
dummy=f1->GetParameter(5);
quadratic=f1->GetParameter(6);
step=f1->GetParameter(7);
TF1 *peak=new TF1("peak",peak,x1,x2,6);
TF1 *bgnd=new TF1("bgnd",background,x1,x2,8);
TF1 *bgndstep=new TF1("bgndstep",backgroundstep,x1,x2,8);
TF1 *bgnd2=new TF1("bgnd2",background2,x1,x2,8);
peak->SetParameters(area,centroid,fwhm,dummy,quadratic,step);
bgnd->SetParameters(bgrstart,bgrslope,area,centroid,fwhm,dummy,quadratic,step);
bgndstep->SetParameters(bgrstart,bgrslope,area,centroid,fwhm,dummy,quadratic,step);
bgnd2->SetParameters(bgrstart,bgrslope,area,centroid,fwhm,dummy,quadratic,step);
peak->SetLineColor(2);
peak->Draw("same");
if(lin==2){
bgnd2->SetLineColor(2);
bgnd2->Draw("same");
}
else{
if(lin==1){
bgnd->SetLineColor(2);
bgnd->Draw("same");
}
if(lin==3){
bgndstep->SetLineColor(2);
bgndstep->Draw("same");
}
}
j=0;
bgr=0;
if(lin==2){
for(i=centroid-fwhm;ipar[3]){
return par[0]+par[1]*(x[0]-par[5]);
}
else{
return par[7]+par[0]+par[1]*(x[0]-par[5]);
}
}
double background(double *x, double *par){
return par[0]+par[1]*(x[0]-par[5]);
}
double background2(double *x, double *par){
return par[0]+par[1]*(x[0]-par[5])+par[6]*(x[0]-par[5])*(x[0]-par[5]);
}