libCPIXE: cpixe-calib/cpixe-calib.c Source File

00001 
00002 /***************************************************************************
00003     Copyright (C) 2007 by Carlos Pascual-Izarra
00004     carlos.pascual_AT_users.sourceforge.net                                                  *
00005 
00006     This program is free software: you can redistribute it and/or modify
00007     it under the terms of the GNU General Public License as published by
00008     the Free Software Foundation, either version 3 of the License, or
00009     (at your option) any later version.
00010 
00011     This program is distributed in the hope that it will be useful,
00012     but WITHOUT ANY WARRANTY; without even the implied warranty of
00013     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00014     GNU General Public License for more details.
00015 
00016     You should have received a copy of the GNU General Public License
00017     along with this program.  If not, see <http://www.gnu.org/licenses/>.
00018 
00019  ***************************************************************************/
00020 
00021 
00032 #ifdef HAVE_CONFIG_H
00033 #include <config.h>
00034 #endif
00035 
00036 #include <stdio.h>
00037 #include <stdlib.h>
00038 #include <math.h>
00039 #include <string.h>
00040 
00041 #include "../libcpixe/libcpixe.h"
00042 #include "../libcpixe/compilopt.h"
00043 
00044 int main(int argc, char *argv[]);
00045 int cpixecalib(void);
00046 
00047 STATFILENAME inputfilename="input-cal.in";
00048 
00056 int main(int argc, char *argv[])
00057 {
00058   int i;
00059   extern STATFILENAME inputfilename;
00060 
00061   if(argc>1){
00062     for(i=1;i<argc;i++){
00063       strcpy(inputfilename,argv[i]);
00064       cpixecalib();
00065     }
00066   }
00067   else cpixecalib();
00068 
00069   return EXIT_SUCCESS;
00070 }
00071 
00080 int cpixecalib(void)
00081 {
00082   extern STATFILENAME inputfilename;
00083   STATFILENAME fckfilename="";
00084   STATFILENAME abscoeffilename="";
00085   STATFILENAME xenerfilename="";
00086 
00087   EXTRAINFO ExtraInfo;
00088   LYR *lyr;
00089   CalibYld *XYldSums;
00090   int i,j,l,Z;
00091   EXP_PARAM *pExpPar;
00092   SIM_PARAM *pSimPar;
00093   CalibYld tempRawEff,tempEff;
00094 
00095    /*Static variables (Those that should be maintained between calls of CPIXEMAIN()*/
00096   static int *PresentElems;
00097   static CPIXERESULTS *CalcFlags;
00098   static CalibYld *LinesEnerArray;
00099   static CalibYld *EffArray;
00100   static FluorCKCoef *FCKCoefArray;
00101   static AbsCoef *TotAbsCoefArray;
00102   static SPT *SPTArray;
00103   static FILTER Filter;
00104   static FILTER DetectorFilter;
00105 
00106   static EXP_PARAM ExpPar;
00107   static foil *Sample;
00108   static int NFoil, NFoilUsed, dimSums;
00109   static FILE *outputfile;
00110 
00111   #if CPIXEVERBOSITY > 0
00112   printf("\n\n\n************ BEGIN OF CPIXE OUTPUT ****************\n\n");
00113   #endif
00114 
00115   pExpPar=&ExpPar;
00116   pSimPar=&pExpPar->simpar;
00117 
00118 
00119   /**********BEGIN INITIALIZATION*********/
00120 
00121   readINPUT(inputfilename, pExpPar, &ExtraInfo);
00122   //Check that a calibration is required
00123 
00124   //For calibration, areas with errors are needed:
00125   if(ExtraInfo.AreasFormat!=2){fprintf(stderr,"\nERROR: Areas + errors are needed for calibration (use DTAREASFILE instead of CALCFLAGS)\n");exit(1);}
00126 
00127   readsample(ExtraInfo.SampleFileNm, &pSimPar->MaxZinsample, &NFoil, &Sample);
00128   readFilter(ExtraInfo.FilterFileNm, &Filter);
00129   readFilter(ExtraInfo.DetectorFileNm, &DetectorFilter);
00130   //translate Filter definition
00131   Filter.changes=1;    //First time at least Filter should be read:
00132   DetectorFilter.changes=1;    //First time at least Filter should be read:
00133 
00134   snprintf(fckfilename, FILENMLENGTH,"%sSC_PRA74.fck",ExtraInfo.DBpath);
00135   snprintf(abscoeffilename, FILENMLENGTH,"%sxabs.abs",ExtraInfo.DBpath);
00136   snprintf(xenerfilename, FILENMLENGTH,"%sxener.dat",ExtraInfo.DBpath);
00137 
00138   #if CPIXEVERBOSITY > 0
00139   printf("\nInitializing. Reading databases...\n\n");
00140   printf("\nStopping Coefs DB: '%sSCOEF.95[A/B]'",ExtraInfo.DBpath);
00141   printf("\nAbsorption Coefs DB: '%s'",abscoeffilename);
00142   printf("\nFluorescence Coefs DB: '%s'",fckfilename);
00143   printf("\nX-ray Energies: '%s'",xenerfilename);
00144   printf("\nDetector Calibration DB: '%s'",ExtraInfo.CalibFileNm);
00145   printf("\nLines to calculate defined in: '%s'\n",ExtraInfo.AreasFileNm);
00146   #endif
00147 
00148   createPresentElems(pSimPar->MaxZinsample, NFoil, Sample, &PresentElems);
00149   readCalcFlags(ExtraInfo.AreasFileNm, PresentElems, pSimPar->MaxZinsample, ExtraInfo.AreasFormat, &CalcFlags);
00150 
00151   //Read the energies for the lines used
00152   readLinesEner(xenerfilename, PresentElems, pSimPar->MaxZinsample, &LinesEnerArray);
00153 
00154   //Instead of reading a regular calib file, read a "fake one" containing all 1's
00155   snprintf(ExtraInfo.CalibFileNm,FILENMLENGTH,"%sidealdet.eff",ExtraInfo.DBpath);
00156 //  readEff(ExtraInfo.CalibFileNm, PresentElems, pSimPar->MaxZinsample, &EffArray);
00157   readEff3(ExtraInfo.CalibFileNm, PresentElems, pSimPar->MaxZinsample, LinesEnerArray, &EffArray); // changed so that cpixecalib uses readEff3 instead of readEff()
00158 
00159 
00160   //Read the fluorescence coefs.
00161   readFCK(fckfilename, pSimPar->MaxZinsample,&FCKCoefArray);
00162 //   readAbsCoef_OLD(abscoeffilename, pSimPar->MaxZinsample, PresentElems, &Filter, &TotAbsCoefArray);
00163   readAbsCoef(abscoeffilename, pSimPar->MaxZinsample, &TotAbsCoefArray);
00164   createSPTs(ExtraInfo.DBpath,pExpPar, pSimPar->MaxZinsample, PresentElems, (double)(2*pExpPar->ion.A), &SPTArray);
00165 
00166   /********** END INITIALIZATION *********/
00167 
00168 
00169 
00170   /******* BEGIN CALCULATION*********/
00171 
00172   //Calculate Filter Transmission if the filter changed
00173   if(Filter.changes){
00174     FilterTrans(pExpPar->simpar.MaxZinsample, LinesEnerArray, TotAbsCoefArray, PresentElems, &Filter);
00176     Filter.changes=0;
00177   }
00178   //Calculate DetectorFilter Transmission if the DetectorFilter changed
00179   if(DetectorFilter.changes){
00180     FilterTrans(pExpPar->simpar.MaxZinsample, LinesEnerArray, TotAbsCoefArray, PresentElems, &DetectorFilter);
00182     DetectorFilter.changes=0;
00183   }
00184 
00185 
00186   initlyrarray(pExpPar, Sample, LinesEnerArray, TotAbsCoefArray, SPTArray, PresentElems, NFoil, &NFoilUsed, &lyr);
00187   dimSums=pSimPar->MaxZinsample+1;
00188   XYldSums=(CalibYld*)calloc(dimSums,sizeof(CalibYld));  //Note that XYldSums is 0-initialized
00189   //integrate_Simpson(pExpPar, TotAbsCoefArray, FCKCoefArray, XYldArray ,NFoilUsed, &Filter, lyr, XYldSums);
00190   integrate_Simpson2(pExpPar, TotAbsCoefArray, FCKCoefArray, EffArray ,NFoilUsed, &Filter, lyr, XYldSums);
00191 
00192   #if CPIXEVERBOSITY > 1
00193   for(Z=1;Z<pSimPar->MaxZinsample+1;Z++){
00194     if(PresentElems[Z]){
00195     fprintf(stdout,"\nXYldResults for Z=%d\n",Z);
00196     fprintCALIBYLD(stdout,&XYldSums[Z]);
00197     }
00198   }
00199   #endif
00200   #if CPIXEVERBOSITY > 1
00201   for(Z=1;Z<pSimPar->MaxZinsample+1;Z++){
00202     if(PresentElems[Z]){
00203     fprintf(stdout,"\nDetector Transmissions for Z=%d\n",Z);
00204     fprintCALIBYLD(stdout,&DetectorFilter.Trans[Z]);
00205     }
00206   }
00207   #endif
00208   #if CPIXEVERBOSITY > 1
00209   for(Z=1;Z<pSimPar->MaxZinsample+1;Z++){
00210     if(PresentElems[Z]){
00211     fprintf(stdout,"\nFilter Transmissions for Z=%d\n",Z);
00212     fprintCALIBYLD(stdout,&Filter.Trans[Z]);
00213     }
00214   }
00215   #endif
00216 
00219     /******END CALCULATION**********/
00220 
00221   /**********BEGIN FILE OUTPUT LINES*******/
00222   if(ExtraInfo.WantOutputfile){
00223     outputfile = fopen(ExtraInfo.OutputFileNm, "w");
00224     if(outputfile== NULL)
00225       {fprintf(stderr,"\nERROR: Could not open file '%s' for write.\n",ExtraInfo.OutputFileNm);exit(1);}
00226 
00227     fprintf(outputfile,"\n\nBEGIN_HEADER \n");
00228     //Write initialization parameters to output file
00229 
00230     fprintf(outputfile,"\nCPIXE-CALIB Output File\n\n");
00231     fprintf(outputfile,"\nStopping Coefs DB: '%sSCOEF.95[A/B]'",ExtraInfo.DBpath);
00232     fprintf(outputfile,"\nAbsorption Coefs DB: '%s'",abscoeffilename);
00233     fprintf(outputfile,"\nFluorescence Coefs DB: '%s'",fckfilename);
00234     fprintf(outputfile,"\nFake Detector Calibration DB: '%s'",ExtraInfo.CalibFileNm);
00235     fprintf(outputfile,"\nExperimental Areas read from: '%s'\n",ExtraInfo.AreasFileNm);
00236 
00237 
00238     fprintf(outputfile,"\n\nKEV      \t%lg\nINCANG   \t%lg\nEXITANG  \t%lg\n",
00239           pExpPar->BeamEner, pExpPar->IncAng/DEG2RAD, pExpPar->DetAng/DEG2RAD);
00240     fprintf(outputfile,"\nCHARGE   \t%lg\nDTCC     \t%lg\n",
00241           pExpPar->simpar.ColCharge, pExpPar->simpar.DTCC);
00242     fprintf(outputfile,"\n\nSANGFRACT    \t%le\nFLUENCE     \t%le\n", pExpPar->SAngFract,pExpPar->Fluence);
00243 
00244     for(i=0;i<NFoil;i++){
00245       fprintf(outputfile,"\nSample[%d]  (%d elem)\n",i,Sample[i].nfoilelm);
00246       fprintf(outputfile,"\n---Eff. thickness=%.3le  (%d sublyrs)\n",lyr[i+1].ThickIn,lyr[i+1].FESxlen);
00247 
00248       for(j=0;j<Sample[i].comp.nelem;j++){
00249       fprintf(outputfile,"'%s'(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
00250           ChemicalSymbol[Sample[i].comp.elem[j].Z],Sample[i].comp.elem[j].Z,j,Sample[i].comp.X[j],100*Sample[i].comp.xn[j],100*Sample[i].comp.w[j]);
00251       }
00252       fprintf(outputfile,"\nFinal Energy in this layer=%lf keV\n\n",lyr[i+1].FoilOutEner);
00253     }
00254 
00255     for(i=0;i<Filter.nlyr;i++){
00256       fprintf(outputfile,"\nFilter[%d] (%lf 1e15at/cm2)  (%d elem)  \n",i,Filter.foil[i].thick,Filter.foil[i].nfoilelm);
00257       for(j=0;j<Filter.foil[i].comp.nelem;j++){
00258       fprintf(outputfile,"(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
00259           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]);
00260       }
00261     }
00262 
00263     fprintf(outputfile,
00264     "Datablocks Format: \n\
00265     Atnumb\n\
00266     \t{emission energies}\n\
00267     \t{Raw Efficiencies}\n\
00268     \t{Corrected Efficiencies}\n");
00269 
00270     fprintf(outputfile,"\n\nEND_HEADER \n");
00271     fprintf(outputfile,"\n\nDATA \n");
00272 
00273     for(Z=1;Z<pSimPar->MaxZinsample+1;Z++){
00274       if(PresentElems[Z]){
00275         //Z2mass(Z, &tempM, 'n');
00276         //Xprod(&TotAbsCoefArray[Z], &FCKCoefArray[Z], pExpPar->ion.Z, Z, tempM , pExpPar->simpar.CalEner, &tempXprodEcal);
00277 
00278         //calculate Calibration Yields (referred to Ecal), Detector efficiencies and Line ratios
00279         for (j = 1; j <= 3; j++){
00280           if (XYldSums[Z].K_[j]>0){
00282             //tempCYld.K_[j]=CalcFlags[Z].simareas.K_[j]/XYldSums[Z].K_[j];
00283             tempRawEff.K_[j]=CalcFlags[Z].simareas.K_[j]/XYldSums[Z].K_[j];
00284             tempEff.K_[j]=tempRawEff.K_[j]/DetectorFilter.Trans[Z].K_[j];
00285             //tempRatio.K_[j]=tempCYld.K_[j]/tempCYld.K_[1];
00286           }
00287           else{
00288             tempRawEff.K_[j]=0;
00289             tempEff.K_[j]=0;
00290             //tempRatio.K_[j]=0;
00291           }
00292         }
00293         for (j = 1; j <= 3; j++) {
00294           for (l = 1; l <= 3; l++){
00295             if (XYldSums[Z].L_[j][l]>0){
00296               //tempCYld.L_[j][l]=CalcFlags[Z].simareas.L_[j][l]/XYldSums[Z].L_[j][l];
00297               tempRawEff.L_[j][l]=CalcFlags[Z].simareas.L_[j][l]/XYldSums[Z].L_[j][l];
00298               tempEff.L_[j][l]=tempRawEff.L_[j][l]/DetectorFilter.Trans[Z].L_[j][l];
00299               //tempRatio.L_[j][l]=tempCYld.L_[j][l]/tempCYld.L_[j][1];
00300             }
00301             else{
00302               tempRawEff.L_[j][l]=0;
00303               tempEff.L_[j][l]=0;
00304               //tempRatio.L_[j][l]=0;
00305             }
00306           }
00307         }
00308         for (j = 1; j <= 3; j++) {
00309           if (XYldSums[Z].M_[j]>0){
00310             //tempCYld.M_[j]=CalcFlags[Z].simareas.M_[j]/XYldSums[Z].M_[j];
00311             tempRawEff.M_[j]=CalcFlags[Z].simareas.M_[j]/XYldSums[Z].M_[j];
00312             tempEff.M_[j]=tempRawEff.M_[j]/DetectorFilter.Trans[Z].M_[j];
00313             //tempRatio.M_[j]=tempCYld.M_[j]/tempCYld.M_[1];
00314           }
00315           else {
00316             tempRawEff.M_[j]=0;
00317             tempEff.M_[j]=0;
00318             //tempRatio.M_[j]=0;
00319           }
00320         }
00321         //Print results:
00322         fprintf(outputfile,"\n %d\n",CalcFlags[Z].atnum); //atomic number
00323         fprintCALIBYLD(outputfile,&LinesEnerArray[Z]);    //emission energies
00324         fprintf(outputfile,"\n");
00325         fprintCALIBYLD(outputfile,&tempRawEff);           //Raw Efficiencies
00326         fprintf(outputfile,"\n");
00327         fprintCALIBYLD(outputfile,&tempEff);              //Corrected efficiencies
00328         fprintf(outputfile,"\n");
00329 //         fprintCALIBYLD(outputfile,&tempRatio);            //Yield ratios (referred to the main line)
00330 //         fprintf(outputfile,"\n");
00331 
00332 
00333         //Calculate and print
00334       }
00335     }
00336   }
00337   /**********END FILE OUTPUT LINES*******/
00338 
00339 
00340   #if CPIXEVERBOSITY > 0
00341   /**********BEGIN DEBUG LINES*******/
00342   /* Just write stuff for debugging*/
00343   printf("\n\n Ebeam=%lf IncAng=%lf  DetAng=%lf \n",
00344         pExpPar->BeamEner, pExpPar->IncAng/DEG2RAD, pExpPar->DetAng/DEG2RAD);
00345   printf("\n Charge=%lf   DTCC=%lf DetColFac=%lf\n",
00346         pExpPar->simpar.ColCharge, pExpPar->simpar.DTCC, pExpPar->DetColFac);
00347   printf("\n CalibEner=%lf\n", pExpPar->simpar.CalEner);
00348 
00349   for(i=0;i<NFoil;i++){
00350     printf("\nSample[%d]  (%d elem)\n",i,Sample[i].nfoilelm);
00351     printf("\n---Eff. thickness=%.3le  (%d sublyrs)\n",lyr[i+1].ThickIn,lyr[i+1].FESxlen);
00352 
00353     for(j=0;j<Sample[i].comp.nelem;j++){
00354     printf("'%s'(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
00355         ChemicalSymbol[Sample[i].comp.elem[j].Z],Sample[i].comp.elem[j].Z,j,Sample[i].comp.X[j],100*Sample[i].comp.xn[j],100*Sample[i].comp.w[j]);
00356     }
00357     printf("\nFinal Energy in this layer=%lf keV\n\n",lyr[i+1].FoilOutEner);
00358   }
00359 
00360   for(i=0;i<Filter.nlyr;i++){
00361     printf("\nFilter[%d] (%lf 1e15at/cm2)  (%d elem)  \n",i,Filter.foil[i].thick,Filter.foil[i].nfoilelm);
00362     for(j=0;j<Filter.foil[i].comp.nelem;j++){
00363     printf("(Z=%d) X[%d]=%lf (%lfat%% / %lfM%%)\n",
00364         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]);
00365     }
00366   }
00367 
00368   for(i=0;i<dimSums;i++){
00369     if(i<=pSimPar->MaxZinsample && PresentElems[i]){
00370       printf("\n\n ************RESULTS for %s (Z=%d)\n",ChemicalSymbol[i],i);
00371       /*K lines:*/
00372       for (j = 1; j <= 3; j++)
00373         if(XYldSums[i].K_[j] > 0. && CalcFlags[i].simareas.K_[j]> 0.)
00374           printf("\n%10s (%2.2lfkeV)\t%le counts (Eff=%le)", LineNm[1][j], LinesEnerArray[i].K_[j],XYldSums[i].K_[j],CalcFlags[i].simareas.K_[j]/XYldSums[i].K_[j]);
00375       /*L-Lines*/
00376       for (j = 1; j <= 3; j++)
00377         for (l = 1; l <= 3; l++)
00378           if(XYldSums[i].L_[j][l] > 0. && CalcFlags[i].simareas.L_[j][l])
00379             printf("\n%10s (%2.2lfkeV)\t%le counts (Eff=%le)", LineNm[j+1][l], LinesEnerArray[i].L_[j][l],XYldSums[i].L_[j][l],CalcFlags[i].simareas.L_[j][l]/XYldSums[i].L_[j][l]);
00380       /*M lines:*/
00381       for (j = 1; j <= 3; j++)
00382         if(XYldSums[i].M_[j] > 0. && CalcFlags[i].simareas.M_[j]> 0.)
00383           printf("\n%10s (%2.2lfkeV)\t%le counts (Eff=%le)",LineNm[5][j], LinesEnerArray[i].M_[j],XYldSums[i].M_[j],CalcFlags[i].simareas.M_[j]/XYldSums[i].M_[j]);
00384     }
00385   }
00386 
00387   /**********END DEBUG LINES*******/
00388   #endif
00389 
00390 
00391   /*It is VERY important to call freeReusable() before starting a new iteration*/
00392   freeReusable(NFoil,&lyr,&Sample,&XYldSums);
00393   //    printf("\n!!!!!!! %d  %d",NCALL,Filter.changes);
00394   freeFilter(&Filter);
00395   freeFilter(&DetectorFilter);
00396 
00397   /******** BEGIN TIDING UP BEFORE FINAL EXIT*********/
00398   /*Free everything for a clean exit*/
00399   #if CPIXEVERBOSITY > 0
00400   printf("\n Cleaning CPIXE-CALIB Memory usage...\n\n");
00401   #endif
00402 
00403   for(i=1; i<=pSimPar->MaxZinsample ;i++){
00404     if(PresentElems[i]){
00405       free(SPTArray[i].S);
00406       free(SPTArray[i].E);
00407       free(SPTArray[i].dSdE);
00408     }
00409   }
00410   free(SPTArray);
00411   free(PresentElems);
00412   free(EffArray);
00413   free(LinesEnerArray);
00414   free(FCKCoefArray);
00415   free(TotAbsCoefArray);
00416   free(CalcFlags);
00417   /******** END TIDING UP BEFORE FINAL EXIT*********/
00418 
00419   if(ExtraInfo.WantOutputfile)fclose(outputfile);
00420   freeExtraInfo(&ExtraInfo);
00421 
00422   #if CPIXEVERBOSITY > 0
00423   printf("\n\n\n************ END OF CPIXE OUTPUT ****************\n\n");
00424   #endif
00425 
00426   return EXIT_SUCCESS;
00427 }
00428 
00429 
00430 
00431 
00432