/*
  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 <http://www.gnu.org/licenses/>.
*/
//
// file: xecon.c
//
//This is a ROOT script that performs analysis of 2D beta-gamma spectra 
//using the net-count method.
//
//
//==============================================================================================================
//Copyright (C) 2008 Anders Ringbom, Swedish Defence Research Agency (FOI)
//E-mail: Anders.Ringbom@foi.se

//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 2 of the License, or 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, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA

// This program uses the data analysis framework ROOT wich is covered by the 
// GNU Lesser General Public License, version 2.1. 
// This license is also supplied in the distribution of this program, together with the complete source
// code for ROOT. 
//===============================================================================================================

//The script xecon requires the following files to be installed in the same directory:
// (you can of course put them wherever you like, but then you have to change xecon accordingly)
//
//nlib.c   -   Routines used by xecon.c
//param.txt  - A parameter file usd by the analysis. 
//             Contains optional ROI definitions, gain calibration 
//             tweak parameters, and the memory effect correction factor.
//<input file> A text file holding the data  directory name and names of the data files to be analysed. 
//
// Two files are produced by xecon.c:
//
// xecon_result.txt      - The result of the analysis. Produced each time the script is excecuted. 
//                         The previous xecon_result.txt is replaced.
// xecon_data.txt        - The main results written on a single line. Results are concatenated to xecon_data.txt
//                         when the script is excecuted. This file can be used to read data into other programs 
//                         for plotting etc. 
//
// Execution of this script: 
// At the ROOT - prompt, write:
//
// .x xecon.c("<input file>",<roi flag =0,1>, <plot flag = 0,1>)
//
//
// <input file> : The name of the input file holding the data directory and the data file names.
// <roi flag>   : =0 use ROIs from original data, =1 Use ROI definitions from param.txt
// <plot flag>  : =0 do not create and plot histograms, =1 create and plot histograms
//
// Example:
//
// .x xecon.c("my_data.txt",0,0)
// This will perform and plot an anlysis using the data files in my_data.txt. Histograms will in this case not 
// be plotted, and ROI definitions from the original data is used.
//
// 
// Example files for param.txt and my_data.txt are icluded in the distribution of this script.
//
// The script starts here:
//
#include "<nlib.c>"

void xecon(char* file, int use_own_roi, int plot){
  //
  // Delete any previously created histograms
  //
   gDirectory->Delete("hsgamma_s");
   gDirectory->Delete("hsbeta_s");
   gDirectory->Delete("hggamma_s");
   gDirectory->Delete("hgbeta_s");
   gDirectory->Delete("hdgamma_s");
   gDirectory->Delete("hdbeta_s");
   gDirectory->Delete("hbg_s");
   gDirectory->Delete("hbg_g");
   gDirectory->Delete("hbg_d");
   gDirectory->Delete("hbgs_gamma_s");
   gDirectory->Delete("hbgs_beta_s");
   gDirectory->Delete("hgamma_g");
   gDirectory->Delete("hbeta_g");
   gDirectory->Delete("hbg_gamma_g");
   gDirectory->Delete("hbg_beta_g");
   gDirectory->Delete("hgamma_d");
   gDirectory->Delete("hbeta_d");
   gDirectory->Delete("hbg_gamma_d");
   gDirectory->Delete("hbg_beta_d"); 
   gDirectory->Delete("hsub_gamma");
   gDirectory->Delete("hsub_beta");
   gDirectory->Delete("hbg_sub");
   gDirectory->Delete("hbg_s_gated");
   gDirectory->Delete("hbg_s_gated2");
   gDirectory->Delete("hbgs_gamma_gated");
   gDirectory->Delete("hbgs_beta_gated");
   gDirectory->Delete("hbgd_gamma_scale_s");
   gDirectory->Delete("hbgd_gamma_scale_g");
   gDirectory->Delete("hbg_g_gated");
   gDirectory->Delete("hbg_g_gated2");
   gDirectory->Delete("hbgg_gamma_gated");
   gDirectory->Delete("hbgg_beta_gated");
   //
   //Version of this script to be written in the result file
   //
   char *version;
   version="xecon 1.2.1 (2010-05-17)"; 
   //
   // Parse filenames from the input file
   //
   char sfile[100],gfile[100],dfile[100];
   GetFilenames(file,sfile,gfile,dfile);
   printf("\n");
   printf("Files read from %s:\n",file);
   printf("\n");
   printf("Sample file: %s\n",sfile);
   printf("Gas bk file: %s\n",gfile);
   printf("Det bk file: %s\n",dfile);  
   //
   //Create and fill 2D beta-gamma histograms 
   //
   int bg_chan_ene_s[4], *bg_hist_s;
   int bg_chan_ene_g[4], *bg_hist_g;
   int bg_chan_ene_d[4], *bg_hist_d;
   bg_hist_s=bg_spectrum(sfile, bg_chan_ene_s);
   TH2F *hbg_s = new TH2F("hbg_s","Beta-gamma", bg_chan_ene_s[1]+1, 0, bg_chan_ene_s[1]+1,bg_chan_ene_s[0]+1, 0, bg_chan_ene_s[0]+1);
   fill2Dhist(bg_hist_s,bg_chan_ene_s,"hbg_s");
   bg_hist_g=bg_spectrum(gfile, bg_chan_ene_g);
   TH2F *hbg_g = new TH2F("hbg_g","Beta-gamma", bg_chan_ene_g[1]+1, 0, bg_chan_ene_g[1]+1,bg_chan_ene_g[0]+1, 0, bg_chan_ene_g[0]+1);
   fill2Dhist(bg_hist_g,bg_chan_ene_g,"hbg_g");
   bg_hist_d=bg_spectrum(dfile, bg_chan_ene_d);
   TH2F *hbg_d = new TH2F("hbg_d","Beta-gamma", bg_chan_ene_d[1]+1, 0, bg_chan_ene_d[1]+1,bg_chan_ene_d[0]+1, 0, bg_chan_ene_d[0]+1);
   fill2Dhist(bg_hist_d,bg_chan_ene_d,"hbg_d");
   //
   // Get calibration tweak parameters from the parameter file param.txt
   //
   double shift[4];
   GetTweak("param.txt",shift);
   printf("\n");
   if (use_own_roi == 0){
   printf("ROI definitions from original data files.\n");
      printf("Tweak parameters: offset gamma: %f tweak gamma: %f ofset beta: %f tweak beta: %f\n",shift[0],shift[1],shift[2],shift[3]);
   }
   else {
   printf("ROI definitions from file param.txt used, tweaking parameters not used.\n");
   }
   printf("\n");
   //  
   //Get number of counts inside ROIs. Use ROI definitions from paramete file if that options is selected, 
   // else use ROI definitons from the original data files. The energy calibration tweak parameters are only applied if ROI
   // definitions from the original data file are used. The original ROI definitions are always used for the detector background file. 
   //
   TCutG *roi_s[30];
   TCutG *roi_g[30];
   TCutG *roi_d[30];
   int roi_counts_sam[30],roi_counts_gas[30],roi_counts_det[30],n_roi;
   if (use_own_roi==1){
   if (plot == 1) {
   n_roi=get_own_roi_counts(sfile,roi_counts_sam, roi_s,"hbg_s",bg_chan_ene_s,1);
   }
   else {
   n_roi=get_own_roi_counts(sfile,roi_counts_sam, roi_s,"hbg_s",bg_chan_ene_s,0);
   }
   get_own_roi_counts(gfile,roi_counts_gas, roi_g,"hbg_g",bg_chan_ene_g,0);
   }
   else{
   if (plot == 1){
     n_roi=get_roi_counts_shift(sfile,roi_counts_sam, roi_s,"hbg_s",bg_chan_ene_s,shift,1);
   }
   else {
    n_roi=get_roi_counts_shift(sfile,roi_counts_sam, roi_s,"hbg_s",bg_chan_ene_s,shift,0);
   }
   get_roi_counts_shift(gfile,roi_counts_gas, roi_g,"hbg_g",bg_chan_ene_g,shift,0);
   }
   get_roi_counts(dfile,roi_counts_det, roi_d,"hbg_d",bg_chan_ene_d,0);
   //
   // Create gated 2D - histograms
   //
    TH2F *hbg_s_gated = new TH2F("hbg_s_gated","Gated Beta-gamma 30 keV", bg_chan_ene_s[1]+1, 0, bg_chan_ene_s[1]+1,bg_chan_ene_s[0]+1, 0, bg_chan_ene_s[0]+1);
    TH2F *hbg_g_gated = new TH2F("hbg_g_gated","Gated Beta-gamma 30 keV", bg_chan_ene_g[1]+1, 0, bg_chan_ene_g[1]+1,bg_chan_ene_g[0]+1, 0, bg_chan_ene_g[0]+1);
    TH2F *hbg_s_gated2 = new TH2F("hbg_s_gated2","Gated Beta-gamma 80 keV", bg_chan_ene_s[1]+1, 0, bg_chan_ene_s[1]+1,bg_chan_ene_s[0]+1, 0, bg_chan_ene_s[0]+1);
    TH2F *hbg_g_gated2 = new TH2F("hbg_g_gated2","Gated Beta-gamma 80 keV ", bg_chan_ene_g[1]+1, 0, bg_chan_ene_g[1]+1,bg_chan_ene_g[0]+1, 0, bg_chan_ene_g[0]+1);
    fill_cut_2Dhist(bg_hist_s,bg_chan_ene_s,"hbg_s_gated" ,roi_s[3]);
    fill_cut_2Dhist(bg_hist_g,bg_chan_ene_g,"hbg_g_gated" ,roi_g[3]);
    fill_cut_2Dhist(bg_hist_s,bg_chan_ene_s,"hbg_s_gated2" ,roi_s[2]);
    fill_cut_2Dhist(bg_hist_g,bg_chan_ene_g,"hbg_g_gated2" ,roi_g[2]);
   //
   // Read acquisition times and dates
   //
   char acq_start_date_s[20],acq_start_time_s[20];    // Data acquisition start date and start time
   char acq_start_date_g[20],acq_start_time_g[20];    // Data acquisition start date and start time
   char acq_start_date_d[20],acq_start_time_d[20];    // Data acquisition start date and start time
   double acq_times[2];                               // Real [0] and live [1] data acquisiton time
   acquisition(sfile,acq_start_date_s, acq_start_time_s,acq_times);
   double sample_livetime = acq_times[1];
   double sample_realtime = acq_times[0];
   int sample_start_time = date_time_conv(acq_start_date_s,acq_start_time_s);
   acquisition(gfile,acq_start_date_g, acq_start_time_g,acq_times);
   double gasbk_livetime = acq_times[1];
   double gasbk_realtime = acq_times[0];
   int gasbk_start_time = date_time_conv(acq_start_date_g,acq_start_time_g);
   acquisition(dfile,acq_start_date_d, acq_start_time_d,acq_times);
   double det_livetime = acq_times[1];
   double det_realtime = acq_times[0];
   //
   //Calculate detector background subtraction factors
   //
   double dsfact =  sample_livetime/det_livetime;
   double dgfact =  gasbk_livetime/det_livetime;
   double gsfact =  sample_livetime/gasbk_livetime;
   //
   //Calculate time difference between sample and gas background
   //
   int tao = sample_start_time - gasbk_start_time;
   //
   //Get xe collection and processing parameters
   //
   char coll_start_date[20],coll_start_time[20];  // Collection start date and start time
   char coll_stop_date[20],coll_stop_time[20];    // Collection stop date and stop time
   double  airvolume;                             // Collected air volume   
   double xe_vol[2];                              // Stable xenon volume [0] with error [1]
   double xe_yield[2];                            // Xenon yield [0] with error [1]
   airvolume=collection(sfile,coll_start_date,coll_start_time,coll_stop_date,coll_stop_time);
   processing(sfile,xe_vol,xe_yield);
   int coll_start_t = date_time_conv(coll_start_date,coll_start_time);
   int coll_stop_t = date_time_conv(coll_stop_date,coll_stop_time);
   int coll_time = coll_stop_t - coll_start_t;
   int proc_time = sample_start_time - coll_stop_t;
   //
   //Calculate gas background subtraction factors
   // 
   double Fac[4];
   Fac[0] = calc_F(Xe133_lam,tao,sample_livetime, sample_realtime,gasbk_livetime,gasbk_realtime);
   Fac[1] = calc_F(Xe135_lam,tao,sample_livetime, sample_realtime,gasbk_livetime,gasbk_realtime);
   Fac[2] = calc_F(Xe131m_lam,tao,sample_livetime, sample_realtime,gasbk_livetime,gasbk_realtime);
   Fac[3] = calc_F(Xe133m_lam,tao,sample_livetime, sample_realtime,gasbk_livetime,gasbk_realtime);
   //
   //Add memory effect correction to gas background subtraction factors
   //
   double memf = GetMemf("param.txt");
   printf("Memory effect correction: %f \n",memf);
   Fac[0] = Fac[0]*memf;
   Fac[1] = Fac[1]*memf;
   Fac[2] = Fac[2]*memf;
   Fac[3] = Fac[3]*memf;
   //
   // Calculate gas background subtraction factor for xe-133. 
   //double sgfact = sample_livetime/gasbk_livetime*Fac[0];
   double sgfact = Fac[0];
   //
   //Get interference terms and efficiencies from the sample data file
   //
   double sample_eff[30],sample_eff_err[30];      // b-g Efficicencies for ROI:s with error, 
   double gasbk_eff[30],gasbk_eff_err[30];      // b-g Efficicencies for ROI:s with error, 
   int roi_eff[30];                         
   double sample_ratio[30],sample_ratio_err[30];                // Subtraction factors between ROIs,
   double gasbk_ratio[30],gasbk_ratio_err[30];
   int roi1_rat[30], roi2_rat[30];                // ROI numbers with subtraction factors in ratio
   bg_efficiency(sfile, roi_eff, sample_eff, sample_eff_err);
   ratios(sfile, roi1_rat, roi2_rat, sample_ratio, sample_ratio_err);
   int nrat = ratios(gfile, roi1_rat, roi2_rat, gasbk_ratio, gasbk_ratio_err);
   //
   //Calcuate xe counts in sample and background
   //
   double sample_net_counts[30],gasbk_net_counts[30];
   double s_xe_count[30],s_xe_var[30],s_xe_lc[30];
   double g_xe_count[30],g_xe_var[30],g_xe_lc[30];
   int xrf_s[2],xrf_g[2];
   get_xe_counts(roi_counts_sam,roi_counts_det,sample_ratio,dsfact,s_xe_count, s_xe_var, s_xe_lc,xrf_s);
   get_xe_counts(roi_counts_gas,roi_counts_det,gasbk_ratio,dgfact,g_xe_count, g_xe_var, g_xe_lc,xrf_g);
   //
   //Calculate net counts in ROI 2-10
   //
   int cas;
   double net_count[30],net_var[30], net_lc[30];
   get_net_counts(s_xe_count,g_xe_count,s_xe_var,g_xe_var,net_count,net_var,net_lc,Fac);
   double xe_conc_roi[30], xe_err_roi[30], kroi[30];
   double xe_act_roi[30], xe_act_err_roi[30], karoi[30];
   //
   // Calculate activity for each ROI
   //
   
   cas = get_act_roi(net_count,net_var,sample_eff,sample_livetime,xe_act_roi,xe_act_err_roi,karoi);
   printf("\n");
   printf("Activity for each ROI\n");
   printf("ROI   \t A (mBq) \t Error (mBq)\n");
   printf("--------------------------------------------------------\n");  
   for (int i=1;i<=10;i++){
     printf("%2d \t %f \t %f\n",i, xe_act_roi[i], xe_act_err_roi[i]);
   }
   printf("--------------------------------------------------------\n"); 
   //
   // Calculate activity concentration for each ROI
   //
   cas = get_conc_roi(net_count,net_var,sample_eff,coll_time,proc_time,sample_livetime,xe_vol,xe_conc_roi,xe_err_roi,kroi);
   printf("\n");
   printf("Concentration for each ROI\n");
   printf("ROI   \t C (mBq/m3) \t Error (mBq/m3)\n");
   printf("--------------------------------------------------------\n");  
   for (int i=2;i<=10;i++){
     printf("%2d \t %f \t %f\n",i, xe_conc_roi[i], xe_err_roi[i]);
   }
   printf("--------------------------------------------------------\n"); 
   //
   //Calculate concentration for each isotope
   //
   double conc[30];
   get_conc(cas, xe_conc_roi,xe_err_roi, conc,coll_time,proc_time,sample_livetime);
   //get_conc(cas, xe_conc_roi,xe_err_roi, s_xe_count, s_xe_lc, g_xe_count, g_xe_lc, conc,coll_time,proc_time,sample_livetime);
   //
   //Caclulate critical limits
   //
   double lc[4];
   get_LC_n(net_count,net_var,kroi,cas,lc,coll_time,proc_time,sample_livetime);
   double lca[5]; //Also for Radon
   get_LCA_n(net_count,net_var,karoi,cas,lca);
   //
   //Calculate MDCs and MDAs
   //
   double mdc[4];
   get_MDC(net_lc, lc, kroi,cas,mdc);
   double mda[5]; //Also fro radon
   get_MDA(net_lc, lca, karoi,cas,mda);
   //
   //write result to standard output
   //
   printf("\n");
   printf("        LCA(mBq) MDA (mBq)\n");
   printf("------------------------------\n");             
   printf("Xe133:  %4.2f       %4.2f\n",lca[0],mda[0]);
   printf("Xe131m: %4.2f       %4.2f\n",lca[1],mda[1]);
   printf("Xe133m: %4.2f       %4.2f\n",lca[2],mda[2]);
   printf("Xe135:  %4.2f       %4.2f\n",lca[3],mda[3]);
   printf("Rn222:  %4.2f       %4.2f\n",lca[4],mda[4]);
   printf("------------------------------\n");
   printf("\n");  
   printf("\n");
   printf("        LC(mBq/m3) MDC (mBq/m3)\n");
   printf("------------------------------\n");             
   printf("Xe133:  %4.2f       %4.2f\n",lc[0],mdc[0]);
   printf("Xe131m: %4.2f       %4.2f\n",lc[1],mdc[1]);
   printf("Xe133m: %4.2f       %4.2f\n",lc[2],mdc[2]);
   printf("Xe135:  %4.2f       %4.2f\n",lc[3],mdc[3]);
   printf("------------------------------\n");
   printf("\n");  
   printf("Isotopic concentrations (mBq/m3):\n");
   printf("------------------------------\n");       
   printf("Xe133:  %4.2f  +/-  %4.2f \n",conc[0],conc[1]);
   printf("Xe131m: %4.2f  +/-  %4.2f \n",conc[2],conc[3]);
   printf("Xe133m: %4.2f  +/-  %4.2f \n",conc[4],conc[5]);
   printf("Xe135:  %4.2f  +/-  %4.2f \n",conc[6],conc[7]);
   printf("------------------------------\n"); 
   //
   //Create result file
   //
   double crs[10],crg[10],crd[10];
   for (i=0;i<10;i++){
     crs[i]= roi_counts_sam[i]/sample_livetime;
     crg[i]= roi_counts_gas[i]/gasbk_livetime;
     crd[i]= roi_counts_det[i]/det_livetime;
   }
   FILE *fout2= fopen("xecon_data.txt","a");
   fprintf(fout2,"\n%s %s\t%4d\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%6.2f\t%4.0f\t%3.0f\t%4.0f\t%3.0f\t%4.2f",coll_start_date,coll_start_time,1,conc[0],conc[1],mdc[0],lc[0],conc[2],conc[3],mdc[1],lc[1],conc[4],conc[5],mdc[2],lc[2],conc[6],conc[7],mdc[3],lc[3],s_xe_count[1],sqrt(s_xe_var[1]),g_xe_count[1],sqrt(g_xe_var[1]),xe_vol[0]);
   fclose(fout2);
   FILE *fout= fopen("xecon_result.txt","w");
   fprintf(fout,"\t\t\tXENON SAMPLE ANALYSIS\n");
   fprintf(fout,"\n");
   fprintf(fout,"\n");
   fprintf(fout,"Created by:%s",version);
   fprintf(fout,"\n");
   fprintf(fout,"Sample:\n");
   fprintf(fout,"%s\n",sfile);
   fprintf(fout,"Gas background:\n");
   fprintf(fout,"%s\n",gfile);
   fprintf(fout,"Detector background:\n");
   fprintf(fout,"%s\n",dfile);
   fprintf(fout,"\n");
   fprintf(fout,"Collection start: %s %s Collection stop: %s %s\n",coll_start_date,coll_start_time,coll_stop_date, coll_stop_time);
   fprintf(fout,"Air volume:       %2.2f m3\n",airvolume);
   fprintf(fout,"Collection time:  %d sec\n",coll_time);
   fprintf(fout,"Processing time:  %d sec\n",proc_time);
   fprintf(fout,"Xenon volume:     %f +/- %f cm3\n",xe_vol[0],xe_vol[1]);
   fprintf(fout,"Sample acq start: %s %s\n", acq_start_date_s, acq_start_time_s);
   fprintf(fout,"Sample acq times: %f (real) %f (live)\n", sample_livetime, sample_realtime);
   fprintf(fout,"Gasbk acq start:  %s %s\n", acq_start_date_g, acq_start_time_g);
   fprintf(fout,"Gasbk acq times:  %f (real) %f (live)\n", gasbk_livetime, gasbk_realtime);
   fprintf(fout,"Detbg acq start:  %s %s\n", acq_start_date_d, acq_start_time_d);
   fprintf(fout,"Detbg acq times:  %f (real) %f (live)\n", det_livetime, det_realtime);
   fprintf(fout,"Det subf  sample: %f\n",dsfact);
   fprintf(fout,"Det subf   gasbk: %f\n",dgfact);
   fprintf(fout,"Time difference between s and g measurements: %d sec\n",tao);
   fprintf(fout,"Gas background subtraction factors:\n");
   fprintf(fout,"F_133: %f F_135: %f F_131m: %f F_133m:  %f\n",Fac[0],Fac[1],Fac[2],Fac[3]);
   fprintf(fout,"Det subf   gasbk: %f\n",dgfact);
   fprintf(fout,"Energy calibration tweaking parameters (offset gamma, tweak gamma, offset  beta, tweak beta): %f, %f, %f, %f\n",shift[0],shift[1],shift[2],shift[3]);
   fprintf(fout,"Memory effect correction: %f\n",memf);
   fprintf(fout,"\n");
   fprintf(fout,"Gross counts in ROI:s:\n");
   if (use_own_roi==1)  { 
    fprintf(fout,"ROI limits are user defined,tweaking parameters are not used. \n");
    }
   fprintf(fout,"ROI   \t Sample \t   Gaskb \t   Detbk\n");
   fprintf(fout,"--------------------------------------------------------\n");  
   for (int i=0;i<n_roi;i++){
     fprintf(fout,"%2d \t %6d \t %6d \t %6d\n",i+1, roi_counts_sam[i],roi_counts_gas[i],roi_counts_det[i]);
   }
   fprintf(fout,"--------------------------------------------------------\n");
   if ((xrf_s[0]==1) && (xrf_s[1]==1)) fprintf(fout,"Radon and Xe-133 found in sample and corrected for in the other ROI:s\n"); 
   if ((xrf_s[0]==0) && (xrf_s[1]==1)) fprintf(fout,"Xe-133 found in sample and corrected for in the other ROI:s\n"); 
   if ((xrf_s[0]==1) && (xrf_s[1]==0)) fprintf(fout,"Radon found in sample and corrected for in the other ROI:s\n");
   if ((xrf_s[0]==0) && (xrf_s[1]==0)) fprintf(fout,"No Radon or 133Xe  found in sample\n"); 
   fprintf(fout,"\n");
   fprintf(fout,"Gross count rates in ROI:s:\n");
   fprintf(fout,"ROI   \t Sample \t   Gaskb \t   Detbk\n");
   fprintf(fout,"--------------------------------------------------------\n");  
   for (int i=0;i<n_roi;i++){
     fprintf(fout,"%2d \t %6f \t %6f \t %6f\n",i+1, crs[i],crg[i],crd[i]);
   }
   fprintf(fout,"--------------------------------------------------------\n");
   fprintf(fout,"\n"); 
   fprintf(fout,"Net counts in sample after corr of det background and interfering isotopes:\n");
   fprintf(fout,"ROI   \t Counts \t   Variance \t   Critical limit\n");
   fprintf(fout,"--------------------------------------------------------\n");  
   for (int i=1;i<=10;i++){
     fprintf(fout,"%2d \t %6d \t %6d \t %6d\n",i, s_xe_count[i], s_xe_var[i],s_xe_lc[i]);
   }
   fprintf(fout,"--------------------------------------------------------\n"); 
   if ((xrf_g[0]==1) && (xrf_g[1]==1)) fprintf(fout,"Radon and Xe-133 found in gasbk and corrected for in the other ROI:s\n"); 
   if ((xrf_g[0]==0) && (xrf_g[1]==1)) fprintf(fout,"Xe-133 found in gasbk and corrected for in the other ROI:s\n"); 
   if ((xrf_g[0]==1) && (xrf_g[1]==0)) fprintf(fout,"Radon found in gasbk and corrected for in the other ROI:s\n");
   if ((xrf_g[0]==0) && (xrf_g[1]==0)) fprintf(fout,"No Radon or 133Xe  found in gas background\n");    
   fprintf(fout,"\n");
   fprintf(fout,"Net counts in gasbg after corr of det background and interfering isotopes:\n");
   fprintf(fout,"ROI   \t Counts \t   Variance \t   Critical limit\n");
   fprintf(fout,"--------------------------------------------------------\n");  
   for (int i=1;i<=10;i++){
     fprintf(fout,"%2d \t %6d \t %6d \t %6d\n",i, g_xe_count[i], g_xe_var[i],g_xe_lc[i]);
   }
   fprintf(fout,"--------------------------------------------------------\n");  

   fprintf(fout,"\n");
   fprintf(fout,"Total net counts in sample after subtraction of gas background:\n");
   fprintf(fout,"ROI   \t Counts \t   Variance \t   Critical limit\n");
   fprintf(fout,"--------------------------------------------------------\n");  
   for (int i=2;i<=10;i++){
     fprintf(fout,"%2d \t %6d \t %6d \t %6d\n",i, net_count[i], net_var[i],net_lc[i]);
   }
   fprintf(fout,"--------------------------------------------------------\n");  
   fprintf(fout,"\n");
  
   fprintf(fout,"\n");
   fprintf(fout,"Concentrations for each ROI\n");
   fprintf(fout,"ROI   \t Conc (mBq/m3) \t Error (mBq/m3)\n");
   fprintf(fout,"--------------------------------------------------------\n");  
   for (int i=1;i<=10;i++){
     fprintf(fout,"%2d \t %f \t %f\n",i, xe_conc_roi[i], xe_err_roi[i]);
   }
   fprintf(fout,"--------------------------------------------------------\n"); 
   if (cas==1){
   fprintf(fout,"Xe-131m and Xe-133m found and corrected for in analysis of Xe-133\n");
   fprintf(fout,"ROI 3,7, and 8 used for calculating the Xe-133 concentration\n");  
   }
   if (cas==2){
   fprintf(fout,"Xe-131m found and corrected for in analysis of Xe-133\n");
   fprintf(fout,"ROI 3,7, and 9 used for calculating the Xe-133 concentration\n");   
   }
   if (cas==3){
   fprintf(fout,"Xe-133m found and corrected for in analysis of Xe-133\n"); 
   fprintf(fout,"ROI 3,8, and 10 used for calculating the Xe-133 concentration\n");  
   }
   if (cas==4){
   fprintf(fout,"No metastables found, no correction peformed for that in the analysis of Xe-133\n");
   fprintf(fout,"ROI 3 and 4 used for calculating the Xe-133 concentration\n");   
   }
   fprintf(fout,"\n");
   fprintf(fout,"Isotopic concentrations (mBq/m3):\n");
   fprintf(fout,"---------------------------------\n");       
   fprintf(fout,"Xe133:  %4.2f  +/-  %4.2f \n",conc[0],conc[1]);
   fprintf(fout,"Xe131m: %4.2f  +/-  %4.2f \n",conc[2],conc[3]);
   fprintf(fout,"Xe133m: %4.2f  +/-  %4.2f \n",conc[4],conc[5]);
   fprintf(fout,"Xe135:  %4.2f  +/-  %4.2f \n",conc[6],conc[7]);
   fprintf(fout,"---------------------------------\n");     
   fprintf(fout,"        LC(mBq/m3) MDC (mBq/m3)\n");
   fprintf(fout,"---------------------------------\n");             
   fprintf(fout,"Xe133:  %4.2f       %4.2f\n",lc[0],mdc[0]);
   fprintf(fout,"Xe131m: %4.2f       %4.2f\n",lc[1],mdc[1]);
   fprintf(fout,"Xe133m: %4.2f       %4.2f\n",lc[2],mdc[2]);
   fprintf(fout,"Xe135:  %4.2f       %4.2f\n",lc[3],mdc[3]);
   fprintf(fout,"---------------------------------\n");       
   fclose(fout);
   if (plot == 1){
   //
   // Create rest of the histograms 
   //
   hbg_s_gated->ProjectionY("hbgs_gamma_gated_30", 1, bg_chan_ene_s[0]);
   hbg_s_gated->ProjectionX("hbgs_beta_gated_30", 1, bg_chan_ene_s[1]);
   hbgs_beta_gated_30->SetTitle("Gated Beta coincidence 30 keV (sample)");
   hbgs_gamma_gated_30->SetTitle("Gated Gamma coincidence 30 keV (sample)");
   hbg_g_gated->ProjectionY("hbgg_gamma_gated_30", 1, bg_chan_ene_g[0]);
   hbg_g_gated->ProjectionX("hbgg_beta_gated_30", 1, bg_chan_ene_g[1]);
   hbgg_beta_gated_30->SetTitle("Gated Beta coincidence 30 keV (gasbk)");
   hbgg_gamma_gated_30->SetTitle("Gated Gamma coincidence 30 keV (gasbk)");
   hbg_s_gated2->ProjectionY("hbgs_gamma_gated_80", 1, bg_chan_ene_s[0]);
   hbg_s_gated2->ProjectionX("hbgs_beta_gated_80", 1, bg_chan_ene_s[1]);
   hbgs_beta_gated_80->SetTitle("Gated Beta coincidence 80 keV (sample)");
   hbgs_gamma_gated_80->SetTitle("Gated Gamma coincidence 80 keV (sample)");
   hbg_g_gated2->ProjectionY("hbgg_gamma_gated_80", 1, bg_chan_ene_g[0]);
   hbg_g_gated2->ProjectionX("hbgg_beta_gated_80", 1, bg_chan_ene_g[1]);
   hbgg_beta_gated_80->SetTitle("Gated Beta coincidence 80 keV (gasbk)");
   hbgg_gamma_gated_80->SetTitle("Gated Gamma coincidence 80 keV (gasbk)");
   hbg_s->ProjectionY("hbgs_gamma", 1, bg_chan_ene_s[0]);
   hbg_s->ProjectionX("hbgs_beta", 1, bg_chan_ene_s[1]);
   hbgs_beta->SetTitle("Beta coincidence");
   hbgs_gamma->SetTitle("Gamma coincidence");
   hbgs_beta->GetXaxis()->SetTitle("Beta energy (channels)");
   hbgs_gamma->GetXaxis()->SetTitle("Gamma energy (channels)");
   hbg_g->ProjectionY("hbgg_gamma", 1, bg_chan_ene_g[0]);
   hbg_g->ProjectionX("hbgg_beta", 1, bg_chan_ene_g[1]);
   hbgg_beta->SetTitle("Beta coincidence");
   hbgg_gamma->SetTitle("Gamma coincidence");
   hbgg_beta->GetXaxis()->SetTitle("Beta energy (channels)");
   hbgg_gamma->GetXaxis()->SetTitle("Gamma energy (channels)");
   hbg_d->ProjectionY("hbgd_gamma", 1, bg_chan_ene_d[0]);
   hbg_d->ProjectionX("hbgd_beta", 1, bg_chan_ene_d[1]);
   hbgd_beta->SetTitle("Beta coincidence");
   hbgd_gamma->SetTitle("Gamma coincidence");
   hbgd_beta->GetXaxis()->SetTitle("Beta energy (channels)");
   hbgd_gamma->GetXaxis()->SetTitle("Gamma energy (channels)");
   TH1F *hsub_gamma = new TH1F("hsub_gamma","Gamma bgr subtracted", bg_chan_ene_s[0]+1, 0, bg_chan_ene_s[0]+1);
   TH1F *hsub_beta = new TH1F("hsub_beta","Beta bgr subtracted", bg_chan_ene_s[0]+1, 0, bg_chan_ene_s[0]+1);
   hsub_gamma->Add(hbgs_gamma,1);
   hsub_beta->Add(hbgs_beta,1);
   hsub_gamma->Add(hbgg_gamma,-sgfact);
   hsub_beta->Add(hbgg_beta,-sgfact);
   TH2F *hbg_sub = new TH2F("hbg_sub","Beta-gamma bgr subtracted", bg_chan_ene_s[1]+1, 0, bg_chan_ene_s[1]+1,bg_chan_ene_s[0]+1, 0, bg_chan_ene_s[0]+1);
   hbg_sub->Add(hbg_s,1);
   TH1F *hbgd_gamma_scale_s = new TH1F("hbgd_gamma_scale_s","", bg_chan_ene_s[0]+1, 0, bg_chan_ene_s[0]+1);
   TH1F *hbgd_gamma_scale_g = new TH1F("hbgd_gamma_scale_g","", bg_chan_ene_s[0]+1, 0, bg_chan_ene_s[0]+1);
   hbgd_gamma_scale_s->Add(hbgd_gamma,dsfact);
   hbgd_gamma_scale_g->Add(hbgd_gamma,dgfact);
   //
   //Create and fill singles histograms
   //
   int g_chan_ene_s[2], *g_hist_s;                    
   int b_chan_ene_s[2], *b_hist_s;     
   int g_chan_ene_g[2], *g_hist_g;                    
   int b_chan_ene_g[2], *b_hist_g;     
   int g_chan_ene_d[2], *g_hist_d;                    
   int b_chan_ene_d[2], *b_hist_d;           
   g_hist_s=g_spectrum(sfile, g_chan_ene_s);
   b_hist_s=b_spectrum(sfile, b_chan_ene_s);
   TH1F *hsgamma_s = new TH1F("hsgamma_s","Gamma singles", g_chan_ene_s[0]+1, 0, g_chan_ene_s[0]+1);
   fill_1D_hist(g_hist_s,g_chan_ene_s,"hsgamma_s");
   TH1F *hsbeta_s = new TH1F("hsbeta_s","Beta singles", b_chan_ene_s[0]+1, 0, b_chan_ene_s[0]+1);
   fill_1D_hist(b_hist_s,b_chan_ene_s,"hsbeta_s");

   g_hist_g=g_spectrum(gfile, g_chan_ene_g);
   b_hist_g=b_spectrum(gfile, b_chan_ene_g);
   TH1F *hggamma_s = new TH1F("hggamma_s","Gamma singles", g_chan_ene_g[0]+1, 0, g_chan_ene_g[0]+1);
   fill_1D_hist(g_hist_g,g_chan_ene_g,"hggamma_s");
   TH1F *hgbeta_s = new TH1F("hgbeta_s","Beta singles", b_chan_ene_g[0]+1, 0, b_chan_ene_g[0]+1);
   fill_1D_hist(b_hist_g,b_chan_ene_g,"hgbeta_s");

   g_hist_d=g_spectrum(dfile, g_chan_ene_d);
   b_hist_d=b_spectrum(dfile, b_chan_ene_d);
   TH1F *hdgamma_s = new TH1F("hdgamma_s","Gamma singles", g_chan_ene_d[0]+1, 0, g_chan_ene_d[0]+1);
   fill_1D_hist(g_hist_d,g_chan_ene_d,"hdgamma_s");
   TH1F *hdbeta_s = new TH1F("hdbeta_s","Beta singles", b_chan_ene_d[0]+1, 0, b_chan_ene_d[0]+1);
   fill_1D_hist(b_hist_d,b_chan_ene_d,"hdbeta_s");
   //
   //Plot histograms and cuts
   //
   // Sample
   //
 char stat[50], dete[50], station[20], detector[20];
 char acq_start_s[50],sample_lt[50];
  TCanvas *Myc1 = new TCanvas("Myc1",sfile,800,800);
  hsgamma_s->SetStats(kFALSE);
  hsbeta_s->SetStats(kFALSE);
  hbg_s->SetStats(kFALSE);
  hbgs_gamma->SetStats(kFALSE);
  hbgs_beta->SetStats(kFALSE);
  Myc1->Divide(2,3);
  Myc1->cd(1);
  TPad *p = gPad;
  p->SetFillColor(41);
  TLatex *text = new TLatex();
  text -> SetTextSize(0.1);
  header(sfile,station, detector);
  sprintf(stat,"Station: %s",station);
  sprintf(dete,"Detector: %s",detector);
  sprintf(acq_start_s,"Acq start: %s %s",acq_start_date_s, acq_start_time_s);
  sprintf(sample_lt,"Measurement livetime (sec): %7.0f",sample_livetime);
  text -> DrawLatex(0.35, 0.85, "Sample");
  text -> SetTextSize(0.06);
  text -> DrawLatex(0.05, 0.7, stat);
  text -> DrawLatex(0.05, 0.63, dete);
  text -> DrawLatex(0.05, 0.56, acq_start_s);
  text -> DrawLatex(0.05, 0.49, sample_lt);
  //fprintf(fout,"Sample acq times: %f (real) %f (live)\n", sample_livetime, sample_realtime);
  Myc1->cd(2);
  hbg_s->Draw("colz");
  for (int i =0;i<n_roi;i++){
    roi_s[i]->Draw();
  }
  Myc1->cd(3);
  hsgamma_s->Draw("");
  Myc1->cd(5);
  hsbeta_s->Draw("");
  Myc1->cd(4);
  hbgs_gamma->Draw("");
  hbgd_gamma_scale_s->SetLineColor(kRed); 
  hbgd_gamma_scale_s->Draw("same");
  Myc1->cd(6);
  hbgs_beta->Draw("");
  //
  // Gas background
  //
 char acq_start_g[50],gasbk_lt[50];
  TCanvas *Myc2 = new TCanvas("Myc2",gfile,800,800);
  hggamma_s->SetStats(kFALSE);
  hgbeta_s->SetStats(kFALSE);
  hbg_g->SetStats(kFALSE);
  hbgg_gamma->SetStats(kFALSE);
  hbgg_beta->SetStats(kFALSE);
  Myc2->Divide(2,3);
  Myc2->cd(1);
  TPad *p = gPad;
  p->SetFillColor(41);
  TLatex *text = new TLatex();
  text -> SetTextSize(0.1);
  text -> DrawLatex(0.25, 0.85, "Gas background");
  sprintf(acq_start_g,"Acq start: %s %s",acq_start_date_g, acq_start_time_g);
  sprintf(gasbk_lt,"Measurement livetime (sec): %7.0f",gasbk_livetime);
  text -> SetTextSize(0.06);
  text -> DrawLatex(0.05, 0.7, stat);
  text -> DrawLatex(0.05, 0.63, dete);
  text -> DrawLatex(0.05, 0.56, acq_start_g);
  text -> DrawLatex(0.05, 0.49, gasbk_lt);
  Myc2->cd(2);
  hbg_g->Draw("colz");
  for (int i =0;i<n_roi;i++){
    roi_g[i]->Draw();
  } 
  Myc2->cd(3);
  hggamma_s->Draw("");
  Myc2->cd(5);
  hgbeta_s->Draw("");
  Myc2->cd(4);
  hbgg_gamma->Draw("");
  hbgd_gamma_scale_g->SetLineColor(kRed); 
  hbgd_gamma_scale_g->Draw("same");
  Myc2->cd(6);
  hbgg_beta->Draw("");
  //
  // Detector background
  //
 char acq_start_d[50],detbk_lt[50];
  TCanvas *Myc3 = new TCanvas("Myc3",dfile,800,800);
  hdgamma_s->SetStats(kFALSE);
  hdbeta_s->SetStats(kFALSE);
  hbg_d->SetStats(kFALSE);
  hbgd_gamma->SetStats(kFALSE);
  hbgd_beta->SetStats(kFALSE);
  Myc3->Divide(2,3);
  Myc3->cd(1);
  TPad *p = gPad;
  p->SetFillColor(41);
  TLatex *text = new TLatex();
  text -> SetTextSize(0.1);
  text -> DrawLatex(0.2, 0.85, "Detector background");
  sprintf(acq_start_d,"Acq start: %s %s",acq_start_date_d, acq_start_time_d);
  sprintf(detbk_lt,"Measurement livetime (sec): %7.0f",det_livetime);
  text -> SetTextSize(0.06);
  text -> DrawLatex(0.05, 0.7, stat);
  text -> DrawLatex(0.05, 0.63, dete);
  text -> DrawLatex(0.05, 0.56, acq_start_d);
  text -> DrawLatex(0.05, 0.49, detbk_lt);
  Myc3->cd(2);
  hbg_d->Draw("colz");
  for (int i =0;i<n_roi-1;i++){
    roi_d[i]->Draw();
  } 
  Myc3->cd(3);
  hdgamma_s->Draw("");
  Myc3->cd(5);
  hdbeta_s->Draw("");
  Myc3->cd(4);
  hbgd_gamma->Draw("");
  Myc3->cd(6);
  hbgd_beta->Draw("");
   //
   // Background subtracted
   //
  char colsta[50],colsto[50],xev[50];
  char xe133s[50],xe131ms[50],xe133ms[50],xe135s[50];
  char mxe133s[50],mxe131ms[50],mxe133ms[50],mxe135s[50];
  sprintf(colsta,"Collection start: %s %s",coll_start_date,coll_start_time);
  sprintf(colsto,"Collection stop: %s %s",coll_stop_date,coll_stop_time);
  text -> DrawLatex(0.05, 0.50, xe133s);   
  sprintf(xev,"Xenon volume: %2.2f ml, yield: %1.3f",xe_vol[0],xe_yield[0]);
  sprintf(xe133s,"Xe133:     %4.2f  +/-  %4.2f",conc[0],conc[1]);
  sprintf(xe131ms,"Xe131m:  %4.2f  +/-  %4.2f",conc[2],conc[3]);
  sprintf(xe133ms,"Xe133m:  %4.2f  +/-  %4.2f",conc[4],conc[5]);
  sprintf(xe135s,"Xe135:     %4.2f  +/-  %4.2f",conc[6],conc[7]);
  sprintf(mxe133s,"Xe133:       %4.2f              %4.2f ",lc[0],mdc[0]);
  sprintf(mxe131ms,"Xe131m:    %4.2f              %4.2f ",lc[1],mdc[1]);
  sprintf(mxe133ms,"Xe133m:    %4.2f              %4.2f ",lc[2],mdc[2]);
  sprintf(mxe135s,"Xe135:       %4.2f              %4.2f ",lc[3],mdc[3]);
      //printf("
  TCanvas *Myc4 = new TCanvas("Myc4","Background subtracted sample",800,800);
  hbg_sub->SetStats(kFALSE);
  hsub_gamma->SetStats(kFALSE);
  hsub_beta->SetStats(kFALSE);
  Myc4->Divide(2,2);
  Myc4->cd(1);
  TPad *p = gPad;
  p->SetFillColor(41);
  TLatex *text = new TLatex();
  //text -> SetTextSize(0.1);
  //text -> DrawLatex(0.35, 0.85, "Sample");
  text -> SetTextSize(0.04);
  text -> DrawLatex(0.05, 0.9, stat);
  text -> DrawLatex(0.05, 0.85, dete);
  text -> DrawLatex(0.05, 0.80, colsta);
  text -> DrawLatex(0.05, 0.75, colsto);
  text -> DrawLatex(0.05, 0.70, xev); 
  text -> DrawLatex(0.05, 0.55, "Xenon concentrations (mBq/m3)");   
  text -> DrawLatex(0.05, 0.50, xe133s); 
  text -> DrawLatex(0.05, 0.45, xe131ms);   
  text -> DrawLatex(0.05, 0.40, xe133ms);   
  text -> DrawLatex(0.05, 0.35, xe135s); 
  text -> DrawLatex(0.05, 0.25, "             LC (mBq/m3)    MDC (mBq/m3)");   
  text -> DrawLatex(0.05, 0.20, mxe133s); 
  text -> DrawLatex(0.05, 0.15, mxe131ms);   
  text -> DrawLatex(0.05, 0.10, mxe133ms);   
  text -> DrawLatex(0.05, 0.05, mxe135s); 
  
  Myc4->cd(2);
  hbg_sub->SetMinimum(0);
  //hbg_sub->SetMaximum(4);
  hbg_sub->Draw("colz");
  for (int i =0;i<n_roi;i++){
    roi_s[i]->Draw();
  }
  Myc4->cd(3);
  hsub_gamma->Draw("");
  Myc4->cd(4);
  hbgs_beta_gated_30->SetStats(kFALSE);
  hbgs_beta_gated_30->Draw("");
  //hsub_beta->Draw("");
   //
   //Gated sample histograms  
   //
  hbg_s_gated->SetStats(kFALSE);
  hbgs_gamma_gated_30->SetStats(kFALSE);
  hbgs_gamma_gated_30->SetStats(kFALSE);
  TCanvas *Myc5 = new TCanvas("Myc5","Gated histograms",600,900);
  Myc5->Divide(2,4);
  Myc5->cd(1);
  hbg_s_gated->Draw("colz");
  for (int i =0;i<n_roi;i++){
    roi_s[i]->Draw();
  }
  Myc5->cd(2);
  hbgs_beta_gated_30->Draw("");
  //hbgg_beta_gated->SetFillColor(kRed);
  //hbgg_beta_gated->Draw("same");
  //
   //Gated gas background histogram 30 keV 
   //
  hbg_g_gated->SetStats(kFALSE);
  hbgg_gamma_gated_30->SetStats(kFALSE);
  hbgg_gamma_gated_30->SetStats(kFALSE);
  Myc5->cd(3);
  hbg_g_gated->Draw("colz");
  for (int i =0;i<n_roi;i++){
    roi_s[i]->Draw();
  }
  Myc5->cd(4);
  hbgg_beta_gated_30->Draw("");
 //
   //Gated sample histogram 80 keV  
   //
  hbg_s_gated2->SetStats(kFALSE);
  hbgs_gamma_gated_80->SetStats(kFALSE);
  hbgs_gamma_gated_80->SetStats(kFALSE);
  Myc5->cd(5);
  hbg_s_gated2->Draw("colz");
  for (int i =0;i<n_roi;i++){
    roi_s[i]->Draw();
  }
  Myc5->cd(6);
  hbgs_beta_gated_80->Draw("");
  //hbgg_beta_gated->SetFillColor(kRed);
  //hbgg_beta_gated->Draw("same");
  //
  //Gated gas background histogram 80 keV 
  //
  hbg_g_gated2->SetStats(kFALSE);
  hbgg_gamma_gated_80->SetStats(kFALSE);
  hbgg_gamma_gated_80->SetStats(kFALSE);
  Myc5->cd(7);
  hbg_g_gated2->Draw("colz");
  for (int i =0;i<n_roi;i++){
    roi_s[i]->Draw();
  }
  Myc5->cd(8);
  hbgg_beta_gated_80->Draw("");
   }
 }