/* 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","<title>", 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","<title>", 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","<title>", 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; i<n; i++){ eg_err[i]= 0; } // Create graphs (with or without errors) if (use_err==0){ TGraph *cal_ec = new TGraph(n,ec,eg); TGraph *cal_ce = new TGraph(n,eg,ec); } else if (use_err==1){ TGraph *cal_ec = new TGraphErrors(n,ec,eg,ec_err,eg_err); TGraph *cal_ce = new TGraphErrors(n,eg,ec,eg_err,ec_err); } else { printf("\n%s\n","Wrong error option in cal_fit (possible: 0 or 1)"); return; } 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)"); 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 <year/month/day OR year-month-day>, char time <hour:min:sec>) // Return value: // int <number of minutes since 2004-01-01> // 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; i<n; i++){ eg_err[i]= 0; } // Create graphs (with or without errors) if (use_err==0){ TGraph *cal_ec = new TGraph(n,ec,eg); TGraph *cal_ce = new TGraph(n,eg,ec); } else if (use_err==1){ TGraph *cal_ec = new TGraphErrors(n,ec,eg,ec_err,eg_err); TGraph *cal_ce = new TGraphErrors(n,eg,ec,eg_err,ec_err); } else { printf("\n%s\n","Wrong error option in cal_fit (possible: 0 or 1)"); return; } // //Perform fits // gStyle->SetOptFit(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;i<n_roi;i++){ // printf("%d %f %f %f %f\n",i+1, roi_b_low[i],roi_b_high[i],roi_g_low[i],roi_g_high[i]); //} //Calculate ROI limits in channels polconv(roi_b_low,roi_b_low_c,pb_ce,n_roi); polconv(roi_b_high,roi_b_high_c,pb_ce,n_roi); polconv(roi_g_low,roi_g_low_c,pg_ce,n_roi); polconv(roi_g_high,roi_g_high_c,pg_ce,n_roi); //printf("ROI limits (channels):\n"); for (int i=0;i<n_roi;i++){ croi_b_low_c[i] = (int) (roi_b_low_c[i]+0.5); croi_b_high_c[i] = (int)(roi_b_high_c[i]+0.5); croi_g_low_c[i] = (int) (roi_g_low_c[i]+0.5); croi_g_high_c[i] = (int)(roi_g_high_c[i]+0.5); //printf("%d %d %d %d %d\n",i+1, croi_b_low_c[i],croi_b_high_c[i],croi_g_low_c[i],croi_g_high_c[i]); } char cutname[60]; for (int i = 0; i < n_roi; i++){ sprintf(cutname,"cut%d%s",i,hname); roi[i] = new TCutG(cutname,5); roi[i]->SetPoint(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;i<n_roi;i++){ ccroi[i]=0; } for (int y=0 ; y < bg_chan_ene[0]; y++) { // loop over gamma for (int x=0 ;x < bg_chan_ene[1]; x++){ // loop over beta cont= hist->GetBinContent(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;i<n_roi;i++){ // printf("%d %f %f %f %f\n",i+1, roi_b_low[i],roi_b_high[i],roi_g_low[i],roi_g_high[i]); //} //Calculate ROI limits in channels polconv(roi_b_low,roi_b_low_c,pb_ce,n_roi); polconv(roi_b_high,roi_b_high_c,pb_ce,n_roi); polconv(roi_g_low,roi_g_low_c,pg_ce,n_roi); polconv(roi_g_high,roi_g_high_c,pg_ce,n_roi); // //Tweak parameters // double offset_g = shift[0]; double tweak_g = shift[1]; double offset_b = shift[2]; double tweak_b = shift[3]; //printf("ROI limits (channels):\n"); for (int i=0;i<n_roi;i++){ //tweak calibration roi_b_low_c[i]= tweak_b*roi_b_low_c[i] + offset_b; roi_b_high_c[i]= tweak_b*roi_b_high_c[i] + offset_b; roi_g_low_c[i]= tweak_g*roi_g_low_c[i] + offset_g; roi_g_high_c[i]= tweak_g*roi_g_high_c[i] + offset_g; // croi_b_low_c[i] = (int) (roi_b_low_c[i]+0.5); croi_b_high_c[i] = (int)(roi_b_high_c[i]+0.5); croi_g_low_c[i] = (int) (roi_g_low_c[i]+0.5); croi_g_high_c[i] = (int)(roi_g_high_c[i]+0.5); //printf("%d %d %d %d %d\n",i+1, croi_b_low_c[i],croi_b_high_c[i],croi_g_low_c[i],croi_g_high_c[i]); } char cutname[30]; for (i = 0; i < n_roi; i++){ sprintf(cutname,"cut%d%s",i,hname); roi[i] = new TCutG(cutname,5); roi[i]->SetPoint(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;i<n_roi;i++){ ccroi[i]=0; } for (int y=0 ; y < bg_chan_ene[0]; y++) { // loop over gamma for (int x=0 ;x < bg_chan_ene[1]; x++){ // loop over beta cont= hist->GetBinContent(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;i<n_roi;i++){ // printf("%d %f %f %f %f\n",i+1, roi_b_low_c[i],roi_b_high_c[i],roi_g_low_c[i],roi_g_high_c[i]); croi_b_low_c[i] = (int) roi_b_low_c[i]; croi_b_high_c[i] = (int) roi_b_high_c[i]; croi_g_low_c[i] = (int) roi_g_low_c[i]; croi_g_high_c[i] = (int) roi_g_high_c[i]; } char cutname[30]; for (int i = 0; i < n_roi; i++){ sprintf(cutname,"cut%d%s",i,hname); roi[i] = new TCutG(cutname,5); roi[i]->SetPoint(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;i<n_roi;i++){ ccroi[i]=0; } for (int y=0 ; y < bg_chan_ene[0]; y++) { // loop over gamma for (int x=0 ;x < bg_chan_ene[1]; x++){ // loop over beta cont= hist->GetBinContent(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;i<centroid+fwhm;i++){ bgr+=bgrstart+j*bgrslope+j*j*quadratic; j++; } } else{ if(lin==1){ for(i=centroid-fwhm;i<centroid+fwhm;i++){ bgr+=bgrstart+j*bgrslope; j++; } } if(lin==3){ for(i=centroid-fwhm;i<centroid;i++){ bgr+=step+bgrstart+j*bgrslope; j++; } for(i=centroid;i<centroid+fwhm;i++){ bgr+=bgrstart+j*bgrslope; j++; } } } printf("\nBackground under peak (centroid +/- fwhm): %.3f\n",bgr); printf("\nPeak-to-background ratio: %.4f\n\n",area/bgr); } double fitflin(double *x, double *par){ return background(x,par) + peak(x,&par[2]); } double fitflinstep(double *x, double *par){ return backgroundstep(x,par) + peak(x,&par[2]); } double fitfquad(double *x, double *par){ return background2(x,par) + peak(x,&par[2]); } double peak(double *x, double *par){ double arg = 0; if (par[2]!=0){ arg=2.35*(x[0]-par[1])/par[2]; } return (1.0667*par[0]/par[2])*TMath::Exp(-0.5*arg*arg); } double backgroundstep(double *x, double *par){ if(x[0]>par[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]); }