IAEAphsp/iaea_phsp/iaea_record.cpp

0001 /*
0002  * Copyright (C) 2006 International Atomic Energy Agency
0003  * -----------------------------------------------------------------------------
0004  *
0005  * Permission is hereby granted, free of charge, to any person obtaining a copy
0006  * of this software and associated documentation files (the "Software"), to deal
0007  * in the Software without restriction, including without limitation the rights
0008  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
0009  * copies of the Software, and to permit persons to whom the Software is furnished
0010  * to do so, subject to the following conditions:
0011  *
0012  * The above copyright notice and this permission notice shall be included in all
0013  * copies or substantial portions of the Software.
0014  *
0015  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
0016  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
0017  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
0018  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
0019  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
0020  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
0021  * THE SOFTWARE.
0022  *
0023  *-----------------------------------------------------------------------------
0024  *
0025  *   AUTHORS:
0026  *
0027  *   Roberto Capote Noy, PhD
0028  *   e-mail: R.CapoteNoy@iaea.org (rcapotenoy@yahoo.com)
0029  *   International Atomic Energy Agency
0030  *   Nuclear Data Section, P.O.Box 100
0031  *   Wagramerstrasse 5, Vienna A-1400, AUSTRIA
0032  *   Phone: +431-260021713; Fax: +431-26007
0033  *
0034  *   Iwan Kawrakow, PhD
0035  *   e-mail iwan@irs.phy.nrc.ca
0036  *   Ionizing Radiation Standards
0037  *   Institute for National Measurement Standards
0038  *   National Research Council of Canada Ottawa, ON, K1A 0R6 Canada
0039  *   Phone: +1-613-993 2197, ext.241; Fax: +1-613-952 9865
0040  *
0041  **********************************************************************************
0042  * For documentation
0043  * see http://www-nds.iaea.org/reports-new/indc-reports/indc-nds/indc-nds-0484.pdf
0044  **********************************************************************************/
0045 //#define DEBUG // Comment to avoid printing for every particle write or read
0046
0047 #if (defined WIN32) || (defined WIN64)
0048 #include <iostream>  // so that namespace std becomes defined
0049 #endif
0050 #include <math.h>
0051 #include <cstdio>
0052
0053 #if !(defined WIN32) && !(defined WIN64)
0054 using namespace std;
0055 #endif
0056
0057 #include "iaea_record.h"
0058
0059 short iaea_record_type::initialize()
0060 {
0061   if(p_file == NULL) {
0062      fprintf(stderr, "\n ERROR: Failed to open Phase Space file \n");
0063      return (FAIL);
0064   }
0065
0066   // Defines i/o logic and variable quantities to be stored
0067   // If the value is zero, then corresponding quantity is fixed
0068   ix = 1;
0069   iy = 1;
0070   iz = 1;
0071   iu = 1;
0072   iv = 1;
0073   iw = 1;
0074   iweight = 1;
0075   // Defines a number of EXTRA long variables stored
0076   // (EGS need 2 for incremental history number and LATCH)
0077   iextralong = 1;
0078   if( iextralong >= NUM_EXTRA_LONG)
0079   {
0080      fprintf(stderr, "\n ERROR: Increase NUM_EXTRA_LONG number in iaea_record.h\n");
0081      return (FAIL);
0082   }
0083   // Defines a number of EXTRA float variables stored (EGS could need 1 for ZLAST)
0084   iextrafloat = 0;
0085   if( iextrafloat >= NUM_EXTRA_FLOAT)
0086   {
0087      fprintf(stderr, "\n ERROR: Increase NUM_EXTRA_FLOAT number in iaea_record.h\n");
0088      return (FAIL);
0089   }
0090
0091   return (OK);
0092 }
0093
0094 short iaea_record_type::write_particle()
0095 {
0096   float floatArray[NUM_EXTRA_FLOAT+7];
0097   IAEA_I32 longArray[NUM_EXTRA_LONG];
0098
0099   char ishort = (char) particle;
0100   if(w < 0) ishort = -ishort; // Sign of w is stored in particle type
0101
0102   if( fwrite(&ishort, sizeof(char),   1, p_file) != 1)
0103   {
0104     fprintf(stderr, "\n ERROR: write_particle: Failed to write particle type\n");
0105     return (FAIL);;
0106   }
0107
0108   int reclength = sizeof(char);
0109
0110   if(IsNewHistory > 0) energy *= (-1); // New history is signaled by negative energy
0111
0112   floatArray[0] = energy;
0113
0114   int i = 0;
0115
0116   if(ix > 0) floatArray[++i] = x;
0117   if(iy > 0) floatArray[++i] = y;
0118   if(iz > 0) floatArray[++i] = z;
0119   if(iu > 0) floatArray[++i] = u;
0120   if(iv > 0) floatArray[++i] = v;
0121   if(iweight > 0) floatArray[++i] = weight;
0122
0123   int j;
0124   for(j=0;j<iextrafloat;j++) floatArray[++i] = extrafloat[j];
0125
0126   reclength += (i+1)*sizeof(float);
0127
0128   if( fwrite(floatArray, sizeof(float), (size_t)(i+1), p_file) != (size_t) (i+1))
0129   {
0130      fprintf(stderr, "\n ERROR: write_particle: Failed to write FLOAT phsp data\n");
0131      return (FAIL);
0132   }
0133
0134   if(iextralong > 0)
0135   {
0136      for(j=0;j<iextralong;j++) longArray[j] = extralong[j];
0137      reclength += iextralong*sizeof(IAEA_I32);
0138      if( fwrite(longArray, sizeof(IAEA_I32), (size_t)iextralong, p_file) != (size_t)iextralong)
0139      {
0140         fprintf(stderr, "\n ERROR: write_particle: Failed to write LONG phsp data\n");
0141         return (FAIL);
0142      }
0143   }
0144
0145   if(reclength == 0) return(FAIL);
0146
0147   #ifdef DEBUG
0148   // charge defined
0149   int iaea_charge[MAX_NUM_PARTICLES]={0,-1,+1,0,+1};
0150   int charge = iaea_charge[particle - 1];
0151
0152   printf("\n Wrote a particle with a record lenght %d",reclength);
0153   printf("\n Q %d E %f X %f Y %f Z %f \n\t u %f v %f w %f W %f Part %d \n",
0154   charge, energy, x, y, z, u, v, w, weight, particle);
0155   if( iextrafloat > 0) printf(" EXTRA FLOATs:");
0156   for(j=0;j<iextrafloat;j++) printf(" F%i %f",j+1,extrafloat[j]);
0157   if( iextralong > 0)  printf(" EXTRA LONGs:");
0158   for(j=0;j<iextralong;j++) printf(" L%i %d", j+1,extralong[j]);
0159   printf("\n");
0160   #endif
0161
0162   return(OK);
0163 }
0164
0165 short iaea_record_type::read_particle()
0166 {
0167   float floatArray[NUM_EXTRA_FLOAT+7];
0168   IAEA_I32 longArray[NUM_EXTRA_LONG+7];
0169   //(MACG) -Wshadow int i,j,is,reclength;
0170   int i,j,l,is,reclength;
0171   char ctmp;
0172
0173   // IAEA_I32 pos = ftell(p_file); // To check file position
0174
0175   if( fread(&ctmp, sizeof(char),   1, p_file) != 1) // particle type is always read
0176   {
0177     fprintf(stderr, "\n ERROR: read_particle: Failed to read particle type\n");
0178     return (FAIL);;
0179   }
0180
0181   particle = (short) ctmp;
0182
0183   is = 1; // getting sign of Z director cosine w
0184   if(particle < 0) {is = -1; particle = -particle;}
0185
0186   reclength = sizeof(char);      // particle type is always read
0187   unsigned int rec_to_read = 1;    // energy is always read
0188
0189   if(ix > 0) rec_to_read++;
0190   if(iy > 0) rec_to_read++;
0191   if(iz > 0) rec_to_read++;
0192   if(iu > 0) rec_to_read++;
0193   if(iv > 0) rec_to_read++;
0194   if(iweight > 0) rec_to_read++;
0195   if(iextrafloat>0) rec_to_read += iextrafloat;
0196
0197   if( fread(floatArray, sizeof(float), rec_to_read, p_file) != rec_to_read)
0198   {
0199     fprintf(stderr, "\n ERROR: read_particle: Failed to read FLOATs \n");
0200     return (FAIL);;
0201   }
0202
0203   reclength += rec_to_read*sizeof(float);
0204
0205
0206   IsNewHistory = 0;
0207   if(floatArray[0]<0) IsNewHistory = 1; // like egsnrc
0208   energy = fabs(floatArray[0]);
0209
0210   i = 0;
0211   if(ix > 0) x = floatArray[++i];
0212   if(iy > 0) y = floatArray[++i];
0213   if(iz > 0) z = floatArray[++i];
0214   if(iu > 0) u = floatArray[++i];
0215   if(iv > 0) v = floatArray[++i];
0216   if(iweight > 0) weight = floatArray[++i];
0217   for(j=0;j<iextrafloat;j++) extrafloat[j] = floatArray[++i];
0218
0219   if(iw > 0)
0220   {
0221       w = 0.f;
0222       double aux = (u*u + v*v);
0223       if (aux<=1.0) w = (float) (is * sqrt((float)(1.0 - aux)));
0224       else
0225       {
0226             aux = sqrt((float)aux);
0227             u /= (float)aux;
0228             v /= (float)aux;
0229       }
0230   }
0231
0232   if(iextralong > 0)
0233   {
0234      if( fread(longArray, sizeof(IAEA_I32), (size_t)iextralong, p_file) != (size_t)iextralong)
0235      {
0236        fprintf(stderr, "\n ERROR: read_particle: Failed to read LONGS\n");
0237        return (FAIL);
0238      }
0239      //(MACG) -Wshadow for(int l=0,j=0;j<iextralong;j++) extralong[j] = longArray[l++];
0240      for(l=0,j=0;j<iextralong;j++) extralong[j] = longArray[l++];
0241      reclength += (iextralong)*sizeof(IAEA_I32);
0242   }
0243
0244   #ifdef DEBUG
0245   // charge defined
0246   int iaea_charge[MAX_NUM_PARTICLES]={0,-1,+1,0,+1};
0247   int charge = iaea_charge[particle - 1];
0248
0249   printf("\n Read a particle with a record lenght %d (New History: %d)",
0250                reclength,IsNewHistory);
0251   printf("\n Q %d E %f X %f Y %f Z %f \n\t u %f v %f w %f W %f Part %d \n",
0252   charge, energy, x, y, z, u, v, w, weight, particle);
0253   if( iextrafloat > 0) printf(" EXTRA FLOATs:");
0254   for(j=0;j<iextrafloat;j++) printf(" F%i %f",j+1,extrafloat[j]);
0255   if( iextralong > 0)  printf(" EXTRA LONGs:");
0256   for(j=0;j<iextralong;j++)  printf(" L%i %d",j+1,extralong[j]);
0257   printf("\n");
0258   #endif
0259   return(reclength);
0260 }