cpixe-calib.c

Go to the documentation of this file.

/***************************************************************************
 *   Copyright (C) 2005 by Carlos Pascual-Izarra                           *
 *   carlos.pascual@itn.pt                                                 *
 *                                                                         *
 *   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     *
 *   (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, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/
 
 
 
/**\file cpixe-calib.c
cpixe-calib is an auxiliary program for creating the detector efficiency calibration files used by programs based on the LibCPIXE library.

cpixe-calib uses also the LibCPIXE library

See the libcpixe.c file for more information.

*/


#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>

#include "libcpixe.h"
#include "compilopt.h"

int main(int argc, char *argv[]);
int cpixecalib(void);

STATFILENAME inputfilename="input-cal.in";

/**Main program for cpixe-calib:
 * The executable accepts a list of command line arguments
 * indicating the input file names
 * If none is specified, ("input-cal.in" is assumed).
 *
 * @return EXIT_SUCCESS if normal finishing
 */
int main(int argc, char *argv[])
{
  int i;
  extern STATFILENAME inputfilename;  
  
  if(argc>1){
    for(i=1;i<argc;i++){
      strcpy(inputfilename,argv[i]);
      cpixecalib();
    }
  }
  else cpixecalib();
  
  return EXIT_SUCCESS;
}

/** Calculates the efficiency coefficient based on experimental data and 
 * (hopefully) accurate description of the sample.
 * The efficiency coefs have the following meaning:
 
 "Number of counts that are expected in the detector for a given X-ray emission line induced by a 2MeV proton beam in an ideally thin target, per unit of irradiation charge (uC), per unit of solid angle (msr) and per 10^15 at/cm2 of the element"
 
 * @return EXIT_SUCCESS if no error ocurred and "-1" if some error happened. 
 */
int cpixecalib(void)
{
  STATFILENAME dbpath="./";
  extern STATFILENAME inputfilename;
  STATFILENAME samplefilenm=""; 
  STATFILENAME filterfilenm=""; 
  STATFILENAME calibfilenm="";
  STATFILENAME areasfilenm="";
  STATFILENAME fckfilename="";
  STATFILENAME abscoeffilename=""; 
  STATFILENAME outputfilenm="";
  
  EXTRAINFO ExtraInfo;
  LYR *lyr;
  CalibYld *XYldSums;
  int i,j,l,Z;
  double tempM;
  EXP_PARAM *pExpPar;
  SIM_PARAM *pSimPar;
  CalibYld tempCYld,tempEff,tempRatio,tempXprodEcal;
  
   /*Static variables (Those that should be maintained between calls of CPIXEMAIN()*/
  static int *PresentElems;
  static CPIXERESULTS *CalcFlags;
  static XrayYield *XYldArray;
  static FluorCKCoef *FCKCoefArray;
  static AbsCoef *TotAbsCoefArray;
  static SPT *SPTArray;
  static FILTER Filter;
  
  static EXP_PARAM ExpPar;
  static foil *Sample;
  static int NFoil, NFoilUsed, dimSums;
  static FILE *outputfile;
  
  #if CPIXEVERBOSITY > 0
  printf("\n\n\n************ BEGIN OF CPIXE OUTPUT ****************\n\n");
  #endif

  pExpPar=&ExpPar;
  pSimPar=&pExpPar->simpar;


  /**********BEGIN INITIALIZATION*********/

  //initialize the file name pointers of the ExtraInfo structure
  ExtraInfo.SampleFileNm=samplefilenm;
  ExtraInfo.FilterFileNm=filterfilenm;
  ExtraInfo.XYldFileNm=calibfilenm;
  ExtraInfo.AreasFileNm=areasfilenm;
  ExtraInfo.DBpath=dbpath;
  ExtraInfo.OutputFileNm=outputfilenm;

  readINPUT(inputfilename, pExpPar, &ExtraInfo);
  //Check that a calibration is required
  
  //For calibration, areas with errors are needed:
  if(ExtraInfo.AreasFormat!=2){fprintf(stderr,"\nERROR: Areas + errors are needed for calibration (use DTAREASFILE instead of CALCFLAGS)\n");exit(1);}
  //also the "DoCalibration Flagg must be on
  if(!ExtraInfo.DoCalibration){fprintf(stderr,"\nERROR: The DoCalibration flag is not 1\n");exit(1);}
  
  readsample(samplefilenm, &pSimPar->MaxZinsample, &NFoil, &Sample);
  readFilter(filterfilenm, &Filter);
  //translate Filter definition     
  Filter.changes=1;    //First time at least Filter should be read:
  
  snprintf(fckfilename, FILENMLENGTH,"%sSC_PRA74.fck",dbpath);
  snprintf(abscoeffilename, FILENMLENGTH,"%sxabs.abs",dbpath); 

  #if CPIXEVERBOSITY > 0  
  printf("\nInitializing. Reading databases...\n\n");
  printf("\nStopping Coefs DB: '%sSCOEF.95[A/B]'",dbpath);
  printf("\nAbsorption Coefs DB: '%s'",abscoeffilename);
  printf("\nFluorescence Coefs DB: '%s'",fckfilename);
  printf("\nDetector Calibration DB: '%s'",calibfilenm);
  printf("\nLines to calculate defined in: '%s'\n",areasfilenm);
  #endif

  createPresentElems(pSimPar->MaxZinsample, NFoil, Sample, &PresentElems);
  readCalcFlags(areasfilenm, PresentElems, pSimPar->MaxZinsample, ExtraInfo.AreasFormat, &CalcFlags);
  
  #if CPIXEVERBOSITY > 0
  for(i=1;i<pSimPar->MaxZinsample+1;i++){
    if(PresentElems[i]){
      printf("\nZ=%d  (%d)\n ",i,CalcFlags[i].atnum);
      for (j = 1; j <= 3; j++) printf("\t%.0lf",CalcFlags[i].simareas.K_[j]);
      for (j = 1; j <= 3; j++) { 
        printf("\n");
        for (l = 1; l <= 3; l++)printf("\t%.0lf",CalcFlags[i].simareas.L_[j][l]);
      }
      printf("\n");
      for (j = 1; j <= 3; j++) printf("\t%.0lf",CalcFlags[i].simareas.M_[j]);
    }
  }
  #endif

  //Instead of reading a regular calib file, read a "fake one" containing all 1's
  snprintf(calibfilenm,FILENMLENGTH,"%sCALIB1.CAL",dbpath);
  readXYld(calibfilenm, PresentElems, CalcFlags, pSimPar->MaxZinsample, &pSimPar->CalEner, &XYldArray);  
  
  readFCK(fckfilename, pSimPar->MaxZinsample,&FCKCoefArray);
  readAbsCoef(abscoeffilename, pSimPar->MaxZinsample, PresentElems, &Filter, &TotAbsCoefArray);
  createSPTs(dbpath,pExpPar, pSimPar->MaxZinsample, PresentElems, (double)(2*pExpPar->ion.A), &SPTArray);
  
  /********** END INITIALIZATION *********/

    

  /******* BEGIN CALCULATION*********/
  
  //Calculate Filter Transmission if the filter changed
  if(Filter.changes){
    FilterTrans(pExpPar, XYldArray, TotAbsCoefArray, PresentElems, &Filter);
    ///@todo This indicates that the filter is not going to be varied (Provisional)
    Filter.changes=0;
  }

  initlyrarray(pExpPar, Sample, XYldArray, TotAbsCoefArray, SPTArray, PresentElems, NFoil, &NFoilUsed, &lyr);
  dimSums=pSimPar->MaxZinsample+1;
  XYldSums=(CalibYld*)calloc(dimSums,sizeof(CalibYld));  //Note that XYldSums is 0-initialized
  integrate_Simpson(pExpPar, TotAbsCoefArray, FCKCoefArray, XYldArray ,NFoilUsed, &Filter, lyr, XYldSums);
  
  
  
  ///@todo Divide the calib yields by Xprod(Ecal) to obtain detector efficiency at the given emission energy. Tabulate [emission_energy,efficiency] --->Output
  ///@todo the detector efficiency can be fitted VS energy. The fit can be interpolated and converted back to calib yields for the elements/lines for which no experimental info was available (i.e. construct the calib file) 
    /******END CALCULATION**********/
  
  /**********BEGIN FILE OUTPUT LINES*******/
  if(ExtraInfo.WantOutputfile){
    outputfile = fopen(outputfilenm, "w");
    if(outputfile== NULL)
      {fprintf(stderr,"\nERROR: Could not open file '%s' for write.\n",outputfilenm);exit(1);}
    
    fprintf(outputfile,"\n\nBEGIN_HEADER \n");
    //Write initialization parameters to output file
    
    fprintf(outputfile,"\nCPIXE-CALIB Output File\n\n");
    fprintf(outputfile,"\nStopping Coefs DB: '%sSCOEF.95[A/B]'",dbpath);
    fprintf(outputfile,"\nAbsorption Coefs DB: '%s'",abscoeffilename);
    fprintf(outputfile,"\nFluorescence Coefs DB: '%s'",fckfilename);
    fprintf(outputfile,"\nFake Detector Calibration DB: '%s'",calibfilenm);
    fprintf(outputfile,"\nExperimental Areas read from: '%s'\n",areasfilenm);
  

    fprintf(outputfile,"\n\nKEV      \t%lg\nINCANG   \t%lg\nEXITANG  \t%lg\n",
          pExpPar->BeamEner, pExpPar->IncAng/DEG2RAD, pExpPar->DetAng/DEG2RAD);
    fprintf(outputfile,"\nCHARGE   \t%lg\nDTCC     \t%lg\nDETCOLFAC \t%lg\n",
          pExpPar->simpar.ColCharge, pExpPar->simpar.DTCC, pExpPar->DetColFac);
    fprintf(outputfile,"\n\nCALENER    \t%lg\n", pExpPar->simpar.CalEner);
    
    for(i=0;i<NFoil;i++){
      fprintf(outputfile,"\nSample[%d]  (%d elem)\n",i,Sample[i].nfoilelm);
      fprintf(outputfile,"\n---Eff. thickness=%.3le  (%d sublyrs)\n",lyr[i+1].ThickIn,lyr[i+1].FESxlen);
      
      for(j=0;j<Sample[i].comp.nelem;j++){
      fprintf(outputfile,"'%s'(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
          XYldArray[Sample[i].comp.elem[j].Z].symb,Sample[i].comp.elem[j].Z,j,Sample[i].comp.X[j],100*Sample[i].comp.xn[j],100*Sample[i].comp.w[j]);
      }
      fprintf(outputfile,"\nFinal Energy in this layer=%lf keV\n\n",lyr[i+1].FoilOutEner);
    }
  
    for(i=0;i<Filter.nlyr;i++){
      fprintf(outputfile,"\nFilter[%d] (%lf 1e15at/cm2)  (%d elem)  \n",i,Filter.foil[i].thick,Filter.foil[i].nfoilelm);
      for(j=0;j<Filter.foil[i].comp.nelem;j++){
      fprintf(outputfile,"(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
          Filter.foil[i].comp.elem[j].Z,j,Filter.foil[i].comp.X[j],100*Filter.foil[i].comp.xn[j],100*Filter.foil[i].comp.w[j]);
      }
    }
    
    fprintf(outputfile,
    "Datablocks Format: \n\
    Atnumb\n\
    \t{emission energies}\n\
    \t{CalibYlds}\n\
    \t{Detector efficiencies}\n\
    \t{Yield Ratios}");

    fprintf(outputfile,"\n\nEND_HEADER \n");
    fprintf(outputfile,"\n\nDATA \n");
    
    for(Z=1;Z<pSimPar->MaxZinsample+1;Z++){
      if(PresentElems[Z]){
        Z2mass(Z, &tempM, 'n');
        Xprod(&TotAbsCoefArray[Z], &FCKCoefArray[Z], pExpPar->ion.Z, Z, tempM , pExpPar->simpar.CalEner, &tempXprodEcal);
        
        //calculate Calibration Yields (referred to Ecal), Detector efficiencies and Line ratios
        for (j = 1; j <= 3; j++)
          if (XYldSums[Z].K_[j]>0){
            tempCYld.K_[j]=CalcFlags[Z].simareas.K_[j]/XYldSums[Z].K_[j]; 
            tempEff.K_[j]=tempCYld.K_[j]/tempXprodEcal.K_[j];
            tempRatio.K_[j]=tempCYld.K_[j]/tempCYld.K_[1];
          }
          else{
           tempCYld.K_[j]=0;
           tempEff.K_[j]=0;
           tempRatio.K_[j]=0;
          }
        for (j = 1; j <= 3; j++) { 
          for (l = 1; l <= 3; l++) 
            if (XYldSums[Z].L_[j][l]>0){
              tempCYld.L_[j][l]=CalcFlags[Z].simareas.L_[j][l]/XYldSums[Z].L_[j][l];
              tempEff.L_[j][l]=tempCYld.L_[j][l]/tempXprodEcal.L_[j][l];
              tempRatio.L_[j][l]=tempCYld.L_[j][l]/tempCYld.L_[j][1];
            }
            else{
           tempCYld.L_[j][l]=0;
           tempEff.L_[j][l]=0;
           tempRatio.L_[j][l]=0;
          }
        }
        for (j = 1; j <= 3; j++) 
          if (XYldSums[Z].M_[j]>0){
            tempCYld.M_[j]=CalcFlags[Z].simareas.M_[j]/XYldSums[Z].M_[j];
            tempEff.M_[j]=tempCYld.M_[j]/tempXprodEcal.M_[j];
            tempRatio.M_[j]=tempCYld.M_[j]/tempCYld.M_[1];  
          }
          else {
           tempCYld.M_[j]=0;
           tempEff.M_[j]=0;
           tempRatio.M_[j]=0;
          }
        
        //Print results:
        fprintf(outputfile,"\n %d\n",CalcFlags[Z].atnum); //atomic number
        fprintCALIBYLD(outputfile,&XYldArray[Z].ener);    //emission energies
        fprintf(outputfile,"\n");       
        fprintCALIBYLD(outputfile,&tempCYld);             //Calibration Ylds (referred to Ecal)
        fprintf(outputfile,"\n");                         
        fprintCALIBYLD(outputfile,&tempEff);              //Detector efficiency
        fprintf(outputfile,"\n");                         
        fprintCALIBYLD(outputfile,&tempRatio);            //Yield ratios (referred to the main line)
        fprintf(outputfile,"\n");                         
        
        
        //Calculate and print 
      }
    }
  }
  /**********END FILE OUTPUT LINES*******/
   
      
  #if CPIXEVERBOSITY > 0  
  /**********BEGIN DEBUG LINES*******/
  /* Just write stuff for debugging*/
  printf("\n\n Ebeam=%lf IncAng=%lf  DetAng=%lf \n",
        pExpPar->BeamEner, pExpPar->IncAng/DEG2RAD, pExpPar->DetAng/DEG2RAD);
  printf("\n Charge=%lf   DTCC=%lf DetColFac=%lf\n",
        pExpPar->simpar.ColCharge, pExpPar->simpar.DTCC, pExpPar->DetColFac);
  printf("\n CalibEner=%lf\n", pExpPar->simpar.CalEner);

  for(i=0;i<NFoil;i++){
    printf("\nSample[%d]  (%d elem)\n",i,Sample[i].nfoilelm);
    printf("\n---Eff. thickness=%.3le  (%d sublyrs)\n",lyr[i+1].ThickIn,lyr[i+1].FESxlen);
    
    for(j=0;j<Sample[i].comp.nelem;j++){
    printf("'%s'(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
        XYldArray[Sample[i].comp.elem[j].Z].symb,Sample[i].comp.elem[j].Z,j,Sample[i].comp.X[j],100*Sample[i].comp.xn[j],100*Sample[i].comp.w[j]);
    }
    printf("\nFinal Energy in this layer=%lf keV\n\n",lyr[i+1].FoilOutEner);
  }

  for(i=0;i<Filter.nlyr;i++){
    printf("\nFilter[%d] (%lf 1e15at/cm2)  (%d elem)  \n",i,Filter.foil[i].thick,Filter.foil[i].nfoilelm);
    for(j=0;j<Filter.foil[i].comp.nelem;j++){
    printf("(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
        Filter.foil[i].comp.elem[j].Z,j,Filter.foil[i].comp.X[j],100*Filter.foil[i].comp.xn[j],100*Filter.foil[i].comp.w[j]);
    }
  }
  
  for(i=0;i<dimSums;i++){
    if(i<=pSimPar->MaxZinsample && PresentElems[i]){
      printf("\n\n ************RESULTS for %s (Z=%d)\n",XYldArray[i].symb,i);
      /*K lines:*/
      for (j = 1; j <= 3; j++) 
        if(XYldSums[i].K_[j] > 0.) 
          printf("\n%10s (%2.2lfkeV)\t%le counts @ %lgkeV   (calib=%le)", LineNm[1][j], XYldArray[i].ener.K_[j],CalcFlags[i].simareas.K_[j]/XYldSums[i].K_[j],pExpPar->simpar.CalEner,XYldArray[i].XYld.K_[j]);
      /*L-Lines*/
      for (j = 1; j <= 3; j++)
        for (l = 1; l <= 3; l++) 
          if(XYldSums[i].L_[j][l] > 0.)
            printf("\n%10s (%2.2lfkeV)\t%le counts @ %lgkeV  (calib=%le)", LineNm[j+1][l], XYldArray[i].ener.L_[j][l], CalcFlags[i].simareas.L_[j][l]/XYldSums[i].L_[j][l],pExpPar->simpar.CalEner,XYldArray[i].XYld.L_[j][l]); 
      /*M lines:*/
      for (j = 1; j <= 3; j++) 
        if(XYldSums[i].M_[j] > 0.) 
          printf("\n%10s (%2.2lfkeV)\t%le counts @ %lgkeV  (calib=%le)",LineNm[5][j], XYldArray[i].ener.M_[j], CalcFlags[i].simareas.M_[j]/XYldSums[i].M_[j],pExpPar->simpar.CalEner,XYldArray[i].XYld.M_[j]);
    }
  }

  /**********END DEBUG LINES*******/
  #endif  
  
   
  /*It is VERY important to call freeReusable() before starting a new iteration*/
  freeReusable(NFoil,&lyr,&Sample,&XYldSums);
  //    printf("\n!!!!!!! %d  %d",NCALL,Filter.changes);
  if(Filter.changes)freeFilter(&Filter);
  
   
  /******** BEGIN TIDING UP BEFORE FINAL EXIT*********/
  /*Free everything for a clean exit*/
  #if CPIXEVERBOSITY > 0  
  printf("\nCleaning CPIXE-CALIB Memory usage...\n\n");
  #endif
  
  //Note: Filter.changes=0 implies that freeFilter was not called in last iteration:
  if(!Filter.changes)freeFilter(&Filter);
  
  for(i=1; i<=pSimPar->MaxZinsample ;i++){
    if(PresentElems[i]){
      free(SPTArray[i].S);
      free(SPTArray[i].E);
      free(SPTArray[i].dSdE);
    }
  }
  free(SPTArray);
  free(PresentElems); 
  free(XYldArray);
  free(FCKCoefArray);
  free(TotAbsCoefArray);
  free(CalcFlags);
  /******** END TIDING UP BEFORE FINAL EXIT*********/
  
  
  if(ExtraInfo.WantOutputfile)fclose(outputfile);

  #if CPIXEVERBOSITY > 0
  printf("\n\n\n************ END OF CPIXE OUTPUT ****************\n\n");
  #endif
  
  return EXIT_SUCCESS;
}

---