Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2026-04-18 07:41:08

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 //
0046 // compile in Linux (copy all sources to a separate directory)
0047 // cc *.cpp -lm -lstdc++ -o test_IAEAphsp
0048 //
0049 // Tested in RED HAT LINUX (gcc,cc,g95,icc) and WINDOWS XP (Microsoft C++ v4.2)
0050 //
0051 #if (defined WIN32) || (defined WIN64)
0052 #include <iostream>  // so that namespace std becomes defined
0053 #endif
0054 #include <cstdio>
0055 #include <cstdlib>
0056 #include <cstring>
0057 #include <cmath>
0058 
0059 #if !(defined WIN32) && !(defined WIN64)
0060 using namespace std;
0061 #endif
0062 
0063 #include "utilities.h"
0064 #include "iaea_header.h"
0065 
0066 int iaea_header_type::read_header ()
0067 {
0068     char line[MAX_STR_LEN];
0069 
0070     if(fheader==NULL)
0071     {
0072       printf("\n ERROR: Unable to open header file \n");
0073         return(FAIL);
0074     }
0075 
0076   // ******************************************************************************
0077   // 1. PHSP format
0078 
0079       /*********************************************/
0080       if ( read_block(line,"FILE_TYPE") == FAIL )
0081       {
0082             printf("\nMandatory keyword FILE_TYPE is not defined in input\n");
0083             return FAIL;
0084       }
0085       else file_type = atoi(line);
0086 
0087       /*********************************************/
0088       if ( read_block(line,"CHECKSUM") == FAIL )
0089       {
0090             printf("\nMandatory keyword CHECKSUM is not defined in input\n");
0091             return FAIL;
0092       }
0093       else {
0094           checksum = (IAEA_I64)atof(line);   //atol(line);  //atol() was not working properly in Windows
0095       }
0096 
0097       /*********************************************/
0098       if ( read_block(line,"RECORD_LENGTH") == FAIL )
0099       {
0100             printf("\nMandatory keyword RECORD_LENGTH is not defined in input\n");
0101             return FAIL;
0102       }
0103       else record_length = atoi(line);
0104 
0105       /*********************************************/
0106       if ( read_block(line,"BYTE_ORDER") == FAIL )
0107       {
0108             printf("\nMandatory keyword BYTE_ORDER is not defined in input\n");
0109             return FAIL;
0110       }
0111       else byte_order = atoi(line);
0112 
0113       /*********************************************/
0114     if( get_blockname(line,"RECORD_CONTENTS") == FAIL)
0115     {
0116        printf("\nMandatory keyword RECORD_CONTENTS is not defined in input\n");
0117          return FAIL;
0118     }
0119 
0120     int i;
0121     for (i=0;i<9;i++) record_contents[i] = 0;
0122 
0123     for (i=0;i<9;i++)
0124     {
0125       if( get_string(fheader,line) == FAIL ) return FAIL;
0126       if( *line == SEGMENT_BEG_TOKEN ) break;
0127       record_contents[i] = atoi(line);
0128     };
0129 
0130     for(i=0;i<record_contents[7];i++)
0131     {
0132       if( get_string(fheader,line) == FAIL ) return FAIL;
0133       if( *line == SEGMENT_BEG_TOKEN ) break;
0134       extrafloat_contents[i] = atoi(line);
0135     }
0136 
0137     for(i=0;i<record_contents[8];i++)
0138     {
0139       if( get_string(fheader,line) == FAIL ) return FAIL;
0140       if( *line == SEGMENT_BEG_TOKEN ) break;
0141       extralong_contents[i] = atoi(line);
0142     }
0143 
0144     /*********************************************/
0145     if( get_blockname(line,"RECORD_CONSTANT") == FAIL)
0146     {
0147        printf("\nMandatory keyword RECORD_CONSTANT is not defined in input\n");
0148          return FAIL;
0149     }
0150 
0151     for (i=0;i<7;i++)
0152     {
0153         record_constant[i] = 32000.f;
0154         if(record_contents[i] > 0) continue;
0155         if( get_string(fheader,line) == FAIL ) return FAIL;
0156         if( *line == SEGMENT_BEG_TOKEN ) break;
0157         record_constant[i] = (float)atof(line);
0158     };
0159 
0160 // ******************************************************************************
0161 // 2. Mandatory description of the phsp
0162 
0163       if ( read_block(coordinate_system_description,"COORDINATE_SYSTEM_DESCRIPTION") == FAIL)
0164       {
0165             printf("\nMandatory keyword COORDINATE_SYSTEM_DESCRIPTION is not defined in input\n");
0166             return FAIL;
0167       }
0168 
0169   if(file_type == 1) // For event generators
0170   {
0171     /*********************************************/
0172       if ( read_block(line,"INPUT_FILE_FOR_EVENT_GENERATOR") == FAIL )
0173       {
0174             printf("\nMandatory keyword INPUT_FILE_FOR_EVENT_GENERATOR is not defined in input\n");
0175             return FAIL;
0176       }
0177   }
0178 
0179   if(file_type == 0) // for phsp files
0180   {
0181       /*********************************************/
0182       if ( read_block(line,"ORIG_HISTORIES") == FAIL )
0183       {
0184             printf("\nMandatory keyword ORIG_HISTORIES is not defined in input\n");
0185             return FAIL;
0186       }
0187       else {
0188           orig_histories = (IAEA_I64)atof(line); //atol(line);  //atol(0 was not working properly in Windows
0189           if( orig_histories == 0) printf(
0190            "\n The number of primary particles (ORIG_HISTORIES) is zero in the HEADER !\n");
0191       }
0192 
0193       /*********************************************/
0194       if ( read_block(line,"PARTICLES") == FAIL )
0195       {
0196             printf("\nMandatory keyword PARTICLES is not defined in input\n");
0197             return FAIL;
0198       }
0199       else nParticles = (IAEA_I64)atof(line);   //atol(line);  //atol(0 was not working properly in Windows
0200 
0201       /*********************************************/
0202       for(int itmp=0;itmp<MAX_NUM_PARTICLES;itmp++) particle_number[itmp]=0;
0203       IAEA_I64 npart;
0204       //atol(line); commented  //atol(0 was not working properly in Windows
0205       if ( read_block(line,"PHOTONS") == OK )
0206             {npart = (IAEA_I64)atof(line); particle_number[0] = npart;}
0207       if ( read_block(line,"ELECTRONS") == OK )
0208             {npart = (IAEA_I64)atof(line); particle_number[1] = npart;}
0209       if ( read_block(line,"POSITRONS") == OK )
0210             {npart = (IAEA_I64)atof(line); particle_number[2] = npart;}
0211       if ( read_block(line,"NEUTRONS") == OK )
0212             {npart = (IAEA_I64)atof(line); particle_number[3] = npart;}
0213       if ( read_block(line,"PROTONS") == OK )
0214             {npart = (IAEA_I64)atof(line); particle_number[4] = npart;}
0215   }
0216 // ******************************************************************************
0217 // 3. Mandatory additional information
0218       if ( read_block(line,"IAEA_INDEX") == FAIL )
0219       {
0220             printf("\nMandatory keyword IAEA_INDEX is not defined in input\n");
0221             return FAIL;
0222       }
0223       else iaea_index = atoi(line);
0224 
0225       /*********************************************/
0226       if ( read_block(title,"TITLE") == FAIL )
0227       {
0228             printf("\nMandatory keyword TITLE is not defined in input\n");
0229             return FAIL;
0230       }
0231 
0232       /*********************************************/
0233       if ( read_block(machine_type,"MACHINE_TYPE") == FAIL)
0234       {
0235             printf("\nMandatory keyword MACHINE_TYPE is not defined in input\n");
0236             return FAIL;
0237       }
0238 
0239       /*********************************************/
0240       if ( read_block(MC_code_and_version,"MONTE_CARLO_CODE_VERSION") == FAIL )
0241       {
0242             printf("\nMandatory keyword MONTE_CARLO_CODE_VERSION is not defined in input\n");
0243             return FAIL;
0244       }
0245 
0246       /*********************************************/
0247       if ( read_block(line,"GLOBAL_PHOTON_ENERGY_CUTOFF") == FAIL )
0248       {
0249             printf("\nMandatory keyword GLOBAL_PHOTON_ENERGY_CUTOFF is not defined in input\n");
0250             return FAIL;
0251       }
0252       else global_photon_energy_cutoff = (float)atof(line);
0253 
0254       /*********************************************/
0255       if ( read_block(line,"GLOBAL_PARTICLE_ENERGY_CUTOFF") == FAIL )
0256       {
0257             printf("\nMandatory keyword GLOBAL_PARTICLE_ENERGY_CUTOFF is not defined in input\n");
0258             return FAIL;
0259       }
0260       else global_particle_energy_cutoff = (float)atof(line);
0261 
0262       /*********************************************/
0263       if ( read_block(transport_parameters,"TRANSPORT_PARAMETERS") == FAIL )
0264       {
0265             printf("\nMandatory keyword TRANSPORT_PARAMETERS is not defined in input\n");
0266             return FAIL;
0267       }
0268 
0269 // ******************************************************************************
0270 // 4. Optional description
0271 
0272       if ( read_block(beam_name,"BEAM_NAME") == FAIL )
0273             printf("ERROR reading BEAM_NAME\n");
0274       if ( read_block(field_size,"FIELD_SIZE") == FAIL )
0275             printf("ERROR reading FIELD_SIZE\n");
0276       if ( read_block(nominal_SSD,"NOMINAL_SSD") == FAIL )
0277             printf("ERROR reading NOMINAL_SSD\n");
0278       if ( read_block(variance_reduction_techniques,
0279                         "VARIANCE_REDUCTION_TECHNIQUES") == FAIL )
0280             printf("VARIANCE_REDUCTION_TECHNIQUES\n");
0281       if ( read_block(initial_source_description,"INITIAL_SOURCE_DESCRIPTION") == FAIL )
0282             printf("INITIAL_SOURCE_DESCRIPTION:\n");
0283 
0284       // Documentation sub-section
0285       /*********************************************/
0286       if ( read_block(MC_input_filename,"MC_INPUT_FILENAME") == FAIL )
0287             printf("MC_INPUT_FILENAME\n");
0288       if ( read_block(published_reference,"PUBLISHED_REFERENCE") == FAIL )
0289             printf("PUBLISHED_REFERENCE\n");
0290       if ( read_block(authors,"AUTHORS") == FAIL ) printf("AUTHORS\n");
0291       if ( read_block(institution,"INSTITUTION") == FAIL ) printf("INSTITUTION\n");
0292       if ( read_block(link_validation,"LINK_VALIDATION") == FAIL )
0293             printf("LINK_VALIDATION\n");
0294       if ( read_block(additional_notes,"ADDITIONAL_NOTES") == FAIL )
0295             printf("ADDITIONAL_NOTES\n");
0296 
0297 // ******************************************************************************
0298 // 5. Statistical information
0299       /*********************************************/
0300     if( get_blockname(line,"STATISTICAL_INFORMATION_PARTICLES") == OK)
0301     {
0302         for(i=0;i<MAX_NUM_PARTICLES;i++)
0303         {
0304               if( get_string(fheader,line) == FAIL ) return FAIL;
0305               if( *line == SEGMENT_BEG_TOKEN ) break;
0306 
0307               if(particle_number[i] == 0) continue;
0308               // -------------------------------------------------------
0309               // The fragment below replaces buggy sscanf() function
0310               int index0 =0, index1 =0, len = strlen(line), icnt =0;
0311               float fbuff[6];
0312               do
0313               {   if(advance(line, &index1, len) != OK) break; // Finding first number
0314                     fbuff[icnt++] = (float)atof(line+index1);
0315               if(advance(line, &index0, len) != OK) break; // Skipping spaces
0316               }   while (icnt < 6);
0317               // -------------------------------------------------------
0318               sumParticleWeight[i] = fbuff[0];
0319               minimumWeight[i]         = fbuff[1];
0320             maximumWeight[i]     = fbuff[2];
0321           averageKineticEnergy[i] = fbuff[3];
0322               minimumKineticEnergy[i] = fbuff[4];
0323               maximumKineticEnergy[i] = fbuff[5];
0324         }
0325     }
0326 
0327     if( get_blockname(line,"STATISTICAL_INFORMATION_GEOMETRY") == OK)
0328       {
0329         minimumX = minimumY = minimumZ = 32000.f;
0330         maximumX = maximumY = maximumZ = 0.f;
0331 
0332         for(i=0;i<3;i++)
0333         {
0334             if(record_contents[i] == 1)
0335             {
0336                   if( get_string(fheader,line) == FAIL ) return FAIL;
0337                   if( *line == SEGMENT_BEG_TOKEN ) break;
0338 
0339                 // -------------------------------------------------------
0340                 // The fragment below replaces buggy sscanf() function
0341                 int index0 =0, index1 =0, len = strlen(line), icnt =0;
0342                 float fbuff[2];
0343                 do
0344                 {   if(advance(line, &index1, len) != OK) break; // Finding first number
0345                       fbuff[icnt++] = (float)atof(line+index1);
0346                 if(advance(line, &index0, len) != OK) break; // Skipping spaces
0347                 } while (icnt < 2);
0348                 // -------------------------------------------------------
0349 
0350                   if( i == 0) {minimumX = fbuff[0]; maximumX  = fbuff[1];}
0351                   if( i == 1) {minimumY = fbuff[0]; maximumY  = fbuff[1];}
0352                   if( i == 2) {minimumZ = fbuff[0]; maximumZ  = fbuff[1];}
0353             }
0354         }
0355       }
0356 
0357     return(OK);
0358 }
0359 
0360 int iaea_header_type::write_blockname(const char *blockname)
0361 {
0362   return
0363   (fprintf(fheader,"%c%s%c\n",SEGMENT_BEG_TOKEN,blockname,SEGMENT_END_TOKEN));
0364 }
0365 
0366 int iaea_header_type::get_blockname(char *line, const char *blockname)
0367 {
0368   char *begptr, *endptr;
0369 
0370   // printf("                       Reading block: %s ...\n",blockname);
0371 
0372   if(fheader==NULL)
0373   {
0374     printf("\n ERROR: Opening header file to Get Block \n"); return(FAIL);
0375   }
0376   rewind(fheader) ;
0377 
0378   while( get_string(fheader,line) == OK )
0379   {
0380     if( *line == SEGMENT_BEG_TOKEN )
0381       {
0382          begptr = line + 1 ;
0383              // printf("                       Read line %s ...\n",line);
0384          endptr = strchr(begptr,SEGMENT_END_TOKEN) ;
0385          if( endptr != NULL ) {
0386            *endptr = '\0' ;
0387                // printf("                       Comparing %s ...\n",begptr);
0388            if( strcmp(begptr,blockname) == 0 ) return(OK) ;
0389          }
0390       }
0391   };
0392   return(FAIL) ;
0393 }
0394 
0395 int iaea_header_type::get_block(char *lineread)
0396 {
0397       int read = FAIL, count = 0;
0398 
0399       char line[MAX_STR_LEN];
0400 
0401       strcpy (lineread,""); // Deleting lineread contents
0402       while( get_string(fheader,line) == OK )
0403     {
0404         if( *line == SEGMENT_BEG_TOKEN ) break;
0405         strcat(lineread+count*MAX_NUMB_LINES,line); count++;
0406         read = OK;
0407     };
0408       return (read);
0409 }
0410 
0411 int iaea_header_type::read_block(char *lineread,const char *blockname)
0412 {
0413     char line[MAX_STR_LEN];
0414     if( get_blockname(line,blockname) != OK) return FAIL;
0415       if( get_block(lineread) != OK) return FAIL;
0416       return OK;
0417 }
0418 
0419 int iaea_header_type::set_record_contents(iaea_record_type *p_iaea_record)
0420 {
0421    int i;
0422    for(i=0;i<9;i++) record_contents[i] = 0;
0423    for(i=0;i<7;i++) record_constant[i] = 32000.f;
0424 
0425    if(p_iaea_record->ix > 0) record_contents[0] = 1;
0426    else record_constant[0] = p_iaea_record->x;
0427 
0428    if(p_iaea_record->iy > 0) record_contents[1] = 1;
0429    else record_constant[1] = p_iaea_record->y;
0430 
0431    if(p_iaea_record->iz > 0) record_contents[2] = 1;
0432    else record_constant[2] = p_iaea_record->z;
0433 
0434    if(p_iaea_record->iu > 0) record_contents[3] = 1;
0435    else record_constant[3] = p_iaea_record->u;
0436 
0437    if(p_iaea_record->iv > 0) record_contents[4] = 1;
0438    else record_constant[4] = p_iaea_record->v;
0439 
0440    if(p_iaea_record->iw > 0) record_contents[5] = 1;
0441    else record_constant[5] = p_iaea_record->w;
0442 
0443    if(p_iaea_record->iweight > 0) record_contents[6] = 1;
0444    else record_constant[6] = p_iaea_record->weight;
0445 
0446    if(p_iaea_record->iextrafloat>0) record_contents[7] = p_iaea_record->iextrafloat;
0447    if(p_iaea_record->iextralong>0) record_contents[8] = p_iaea_record->iextralong;
0448 
0449    record_length = 5; // To consider for particle type (1 byte) and energy (4 bytes)
0450    for(i=0;i<8;i++) record_length += record_contents[i]*sizeof(float);
0451    record_length -= 4; // 4 bytes substracted as w is not stored, just his sign
0452    record_length += record_contents[8]*sizeof(IAEA_I32);
0453    if(record_length > 0) return OK;
0454    else
0455    {
0456          printf("\nRECORD LENGTH IS ZERO, CHECK HEADER BLOCK RECORD_CONTENTS\n");
0457          return FAIL;
0458    }
0459 }
0460 
0461 int iaea_header_type::get_record_contents(iaea_record_type *p_iaea_record)
0462 {
0463    int i;
0464 
0465    p_iaea_record->ix = 0;
0466    if(record_contents[0] == 1) p_iaea_record->ix = 1;
0467    else p_iaea_record->x = record_constant[0];
0468 
0469    p_iaea_record->iy = 0;
0470    if(record_contents[1] == 1) p_iaea_record->iy = 1;
0471    else p_iaea_record->y = record_constant[1];
0472 
0473    p_iaea_record->iz = 0;
0474    if(record_contents[2] == 1) p_iaea_record->iz = 1;
0475    else p_iaea_record->z = record_constant[2];
0476 
0477    p_iaea_record->iu = 0;
0478    if(record_contents[3] == 1) p_iaea_record->iu = 1;
0479    else p_iaea_record->u = record_constant[3];
0480 
0481    p_iaea_record->iv = 0;
0482    if(record_contents[4] == 1) p_iaea_record->iv = 1;
0483    else p_iaea_record->v = record_constant[4];
0484 
0485    p_iaea_record->iw = 0;
0486    if(record_contents[5] == 1) p_iaea_record->iw = 1;
0487    else p_iaea_record->w = record_constant[5];
0488 
0489    p_iaea_record->iweight = 0;
0490    if(record_contents[6] == 1) p_iaea_record->iweight = 1;
0491    else p_iaea_record->weight = record_constant[6];
0492 
0493    p_iaea_record->iextrafloat = 0;
0494    if(record_contents[7] > 0) p_iaea_record->iextrafloat = record_contents[7];
0495 
0496    p_iaea_record->iextralong = 0;
0497    if(record_contents[8] > 0) p_iaea_record->iextralong = record_contents[8];
0498 
0499    record_length = 5; // To consider for particle type (1 bytes) and energy (4 bytes)
0500    for(i=0;i<8;i++) record_length += record_contents[i]*sizeof(float);
0501    record_length -= 4; // 4 bytes substracted as w is not stored, just his sign
0502    record_length += record_contents[8]*sizeof(IAEA_I32);
0503    if(record_length > 0) return OK;
0504    else
0505    {
0506          printf("\nWRONG DEFINED HEADER BLOCK RECORD_CONTENTS\n");
0507          return FAIL;
0508    }
0509 }
0510 
0511 void iaea_header_type::initialize_counters()
0512 {
0513   nParticles = read_indep_histories = 0;
0514 
0515   for(int i=0;i<MAX_NUM_PARTICLES;i++)
0516   {
0517         particle_number[i] = 0;
0518         sumParticleWeight[i] = 0.;
0519         maximumKineticEnergy[i] = 0.;
0520         minimumKineticEnergy[i] = 32000.;
0521         averageKineticEnergy[i] = 0.;
0522         minimumWeight[i] = 32000.;
0523         maximumWeight[i] = 0.;
0524   }
0525   minimumX = minimumY = minimumZ = 32000.f;
0526   maximumX = maximumY = maximumZ = -32000.f;
0527 
0528 }
0529 
0530 void iaea_header_type::update_counters(iaea_record_type *p_iaea_record)
0531 {
0532   if (p_iaea_record->x > maximumX )  maximumX = p_iaea_record->x;
0533   if (p_iaea_record->x < minimumX )  minimumX = p_iaea_record->x;
0534 
0535   if (p_iaea_record->y > maximumY )  maximumY = p_iaea_record->y;
0536   if (p_iaea_record->y < minimumY )  minimumY = p_iaea_record->y;
0537 
0538   if (p_iaea_record->z > maximumZ )  maximumZ = p_iaea_record->z;
0539   if (p_iaea_record->z < minimumZ )  minimumZ = p_iaea_record->z;
0540 
0541 
0542   nParticles++;
0543 
0544   if ( p_iaea_record->IsNewHistory > 0 )
0545       read_indep_histories += p_iaea_record->IsNewHistory;
0546 
0547   int i = p_iaea_record->particle-1;
0548   if( i >= 0 && i < MAX_NUM_PARTICLES ) {
0549       particle_number[i]++;
0550       sumParticleWeight[i] +=  p_iaea_record->weight;
0551       averageKineticEnergy[i] += p_iaea_record->weight*
0552                                  fabs(p_iaea_record->energy);
0553       if (p_iaea_record->weight > maximumWeight[i] )
0554             maximumWeight[i] = p_iaea_record->weight;
0555       if (p_iaea_record->weight < minimumWeight[i] )
0556             minimumWeight[i] = p_iaea_record->weight;
0557 
0558       if (fabs(p_iaea_record->energy) > maximumKineticEnergy[i] )
0559          maximumKineticEnergy[i] = fabs(p_iaea_record->energy);
0560       if (fabs(p_iaea_record->energy) < minimumKineticEnergy[i] )
0561          minimumKineticEnergy[i] = fabs(p_iaea_record->energy);
0562   }
0563 
0564 }
0565 
0566 void iaea_header_type::print_statistics()
0567 {
0568    printf("\n *************************************** \n");
0569    printf("           IAEA PHSP STATISTICS          \n");
0570    printf(" *************************************** \n");
0571 
0572 
0573    printf("\n Number of primary particles: %llu\n",orig_histories);
0574    printf(" Number of statistically independent histories read so far %llu\n",read_indep_histories);
0575    printf(" Total number of particles: %llu\n\n",nParticles);
0576 
0577    if(particle_number[0]>0)
0578    {
0579          printf(" PHOTONS: %llu,",particle_number[0]);
0580          if(sumParticleWeight[0]>0) printf(" <E> = %10.4G,",
0581 
0582              averageKineticEnergy[0]/sumParticleWeight[0]);
0583 
0584          printf("  Min.KE = %10.4G, Max.KE = %10.4G\n",
0585 
0586              minimumKineticEnergy[0],maximumKineticEnergy[0]);
0587 
0588          printf("  Total Weight = %15.6G \n",sumParticleWeight[0]);
0589 
0590          printf("  Min.Weight = %12.6G, Max.Weight = %12.6G\n\n",
0591 
0592              minimumWeight[0],maximumWeight[0]);
0593 
0594    }
0595    if(particle_number[1]>0)
0596    {
0597          printf(" ELECTRONS: %llu,",particle_number[1]);
0598          if(sumParticleWeight[1]>0) printf(" <E> = %10.4G,",
0599 
0600              averageKineticEnergy[1]/sumParticleWeight[1]);
0601 
0602          printf("  Min.KE = %10.4G, Max.KE = %10.4G\n",
0603 
0604              minimumKineticEnergy[1],maximumKineticEnergy[1]);
0605 
0606          printf("  Total Weight = %15.6G \n",sumParticleWeight[1]);
0607 
0608          printf("  Min.Weight = %12.6G, Max.Weight = %12.6G\n\n",
0609 
0610              minimumWeight[1],maximumWeight[1]);
0611 
0612    }
0613    if(particle_number[2]>0)
0614    {
0615          printf(" POSITRONS: %llu,",particle_number[2]);
0616          if(sumParticleWeight[2]>0) printf(" <E> = %10.4G,",
0617 
0618              averageKineticEnergy[2]/sumParticleWeight[2]);
0619 
0620          printf("  Min.KE = %10.4G, Max.KE = %10.4G\n",
0621 
0622              minimumKineticEnergy[2],maximumKineticEnergy[2]);
0623 
0624          printf("  Total Weight = %15.6G \n",sumParticleWeight[2]);
0625 
0626          printf("  Min.Weight = %12.6G, Max.Weight = %12.6G\n\n",
0627 
0628              minimumWeight[2],maximumWeight[2]);
0629 
0630    }
0631    if(particle_number[3]>0)
0632    {
0633          printf(" NEUTRONS: %llu,",particle_number[3]);
0634          if(sumParticleWeight[3]>0) printf(" <E> = %10.4G,",
0635 
0636              averageKineticEnergy[3]/sumParticleWeight[3]);
0637 
0638          printf("  Min.KE = %10.4G, Max.KE = %10.4G\n",
0639 
0640              minimumKineticEnergy[3],maximumKineticEnergy[3]);
0641 
0642          printf("  Total Weight = %15.6G \n",sumParticleWeight[3]);
0643 
0644          printf("  Min.Weight = %12.6G, Max.Weight = %12.6G\n\n",
0645 
0646              minimumWeight[3],maximumWeight[3]);
0647 
0648    }
0649    if(particle_number[4]>0)
0650    {
0651          printf(" PROTONS: %llu,",particle_number[4]);
0652          if(sumParticleWeight[4]>0) printf(" <E> = %10.4G,",
0653 
0654              averageKineticEnergy[4]/sumParticleWeight[4]);
0655 
0656          printf("  Min.KE = %10.4G, Max.KE = %10.4G\n",
0657 
0658              minimumKineticEnergy[4],maximumKineticEnergy[4]);
0659 
0660          printf("  Total Weight = %15.6G \n",sumParticleWeight[4]);
0661 
0662          printf("  Min.Weight = %12.6G, Max.Weight = %12.6G\n\n",
0663 
0664              minimumWeight[4],maximumWeight[4]);
0665 
0666    }
0667    printf(" GEOMETRY STATISTICS\n");
0668    if(record_contents[0] == 1)
0669          printf("  %7.3f < X coordinate < %7.3f\n",minimumX,maximumX);
0670    else printf("  X (constant) = %7.3f",record_constant[0]);
0671    if(record_contents[1] == 1)
0672          printf("  %7.3f < Y coordinate < %7.3f\n",minimumY,maximumY);
0673    else printf("  Y (constant) = %7.3f",record_constant[1]);
0674    if(record_contents[2] == 1)
0675          printf("  %7.3f < Z coordinate < %7.3f\n",minimumZ,maximumZ);
0676    else printf("  Z (constant) = %7.3f",record_constant[2]);
0677    printf("\n *************************************** \n\n");
0678 }
0679 
0680 int iaea_header_type::check_byte_order()
0681 {
0682   /* Determine the byte order on this machine */
0683   float ftest=1.0f; /* assign a float to 1.0 */
0684   char *pf = (char *) &ftest;
0685   // printf("\n \t %x %x %x %x", pf[0],pf[1],pf[2],pf[3]);
0686   if(pf[0] == 0 && pf[3] != 0)
0687   {
0688     // printf("\nByte order: INTEL / ALPHA,LINUX -> LITLE_ENDIAN \n");
0689     return(LITTLE_ENDIAN);
0690   }else if(pf[0] != 0 && pf[3] == 0)
0691   {
0692     // printf("\nByte order: OTHER (SGI,SUN-SOLARIS) -> BIG_ENDIAN \n ");
0693     return(BIG_ENDIAN);
0694   }
0695   else
0696   {
0697     printf("\nERROR: indeterminate byte order");
0698     printf("\n \t %x %x %x %x", pf[0],pf[1],pf[2],pf[3]);
0699     return(UNKNOWN_ENDIAN);
0700   }
0701 }
0702 
0703 int iaea_header_type::write_header()
0704 {
0705   if(fheader==NULL)
0706   {
0707       printf("\n ERROR: Opening header file to write \n"); return(FAIL);
0708   }
0709 
0710   rewind(fheader);
0711 
0712   if( write_blockname("IAEA_INDEX") == FAIL ) return(FAIL);
0713 
0714   fprintf(fheader,"%i   // Test header\n\n",iaea_index);
0715 
0716   write_blockname("TITLE");fprintf(fheader,"%s \n\n",title);
0717 
0718   write_blockname("FILE_TYPE");fprintf(fheader,"0\n\n"); // phasespace is assumed
0719 
0720   checksum = (IAEA_I64)record_length * nParticles;
0721 
0722   write_blockname("CHECKSUM");fprintf(fheader,"%llu \n\n",checksum);
0723 
0724   write_blockname("RECORD_CONTENTS");
0725   fprintf(fheader,"   %2i     // X is stored ?\n",record_contents[0]);
0726   fprintf(fheader,"   %2i     // Y is stored ?\n",record_contents[1]);
0727   fprintf(fheader,"   %2i     // Z is stored ?\n",record_contents[2]);
0728 
0729   fprintf(fheader,"   %2i     // U is stored ?\n",record_contents[3]);
0730   fprintf(fheader,"   %2i     // V is stored ?\n",record_contents[4]);
0731   fprintf(fheader,"   %2i     // W is stored ?\n",record_contents[5]);
0732 
0733   fprintf(fheader,"   %2i     // Weight is stored ?\n",record_contents[6]);
0734   fprintf(fheader,"   %2i     // Extra floats stored ?\n",record_contents[7]);
0735   fprintf(fheader,"   %2i     // Extra longs stored ?\n",record_contents[8]);
0736 
0737   int i;
0738   for(i=0;i<record_contents[7];i++) {
0739       if(extrafloat_contents[i] == 0) fprintf(fheader,
0740         "   %2i     // Generic float variable stored in the extrafloat array [%2i] \n",
0741         extrafloat_contents[i], i);
0742       if(extrafloat_contents[i] == 1) fprintf(fheader,
0743         "   %2i     // XLAST variable stored in the extrafloat array [%2i] \n",
0744         extrafloat_contents[i], i);
0745       if(extrafloat_contents[i] == 2) fprintf(fheader,
0746         "   %2i     // YLAST variable stored in the extrafloat array [%2i] \n",
0747         extrafloat_contents[i], i);
0748       if(extrafloat_contents[i] == 3) fprintf(fheader,
0749         "   %2i     // ZLAST variable stored in the extrafloat array [%2i] \n",
0750         extrafloat_contents[i], i);
0751   }
0752 
0753   for(i=0;i<record_contents[8];i++) {
0754       if(extralong_contents[i] == 0) fprintf(fheader,
0755         "   %2i     // Generic integer variable stored in the extralong array [%2i] \n",
0756         extralong_contents[i], i);
0757 
0758       if(extralong_contents[i] == 1) {
0759         fprintf(fheader,
0760         "   %2i     // Incremental history number stored in the extralong array [%2i] \n",
0761         extralong_contents[i], i);
0762         // orig_histories = read_indep_histories;
0763       }
0764       if(extralong_contents[i] == 2) fprintf(fheader,
0765         "   %2i     // LATCH EGS variable stored in the extralong array [%2i] \n",
0766         extralong_contents[i], i);
0767 
0768       if(extralong_contents[i] == 3) fprintf(fheader,
0769         "   %2i     // ILB5 PENELOPE variable stored in the extralong array [%2i] \n",
0770         extralong_contents[i], i);
0771 
0772       if(extralong_contents[i] == 4) fprintf(fheader,
0773         "   %2i     // ILB4 PENELOPE variable stored in the extralong array [%2i] \n",
0774         extralong_contents[i], i);
0775 
0776       if(extralong_contents[i] == 5) fprintf(fheader,
0777         "   %2i     // ILB3 PENELOPE variable stored in the extralong array [%2i] \n",
0778         extralong_contents[i], i);
0779 
0780       if(extralong_contents[i] == 6) fprintf(fheader,
0781         "   %2i     // ILB2 PENELOPE variable stored in the extralong array [%2i] \n",
0782         extralong_contents[i], i);
0783 
0784       if(extralong_contents[i] == 7) fprintf(fheader,
0785         "   %2i     // ILB1 PENELOPE variable stored in the extralong array [%2i] \n",
0786         extralong_contents[i], i);
0787   }
0788 
0789   fprintf(fheader,"\n");
0790 
0791   write_blockname("RECORD_CONSTANT");
0792   if(record_contents[0]==0)
0793     fprintf(fheader,"   %8.4f     // Constant X\n",record_constant[0]);
0794   if(record_contents[1]==0)
0795     fprintf(fheader,"   %8.4f     // Constant Y\n",record_constant[1]);
0796   if(record_contents[2]==0)
0797     fprintf(fheader,"   %8.4f     // Constant Z\n",record_constant[2]);
0798   if(record_contents[3]==0)
0799     fprintf(fheader,"   %8.5f     // Constant U\n",record_constant[3]);
0800   if(record_contents[4]==0)
0801     fprintf(fheader,"   %8.5f     // Constant V\n",record_constant[4]);
0802   if(record_contents[5]==0)
0803     fprintf(fheader,"   %8.5f     // Constant W\n",record_constant[5]);
0804   if(record_contents[6]==0)
0805     fprintf(fheader,"   %8.4f     // Constant Weight\n",record_constant[6]);
0806 
0807   fprintf(fheader,"\n");
0808 
0809   write_blockname("RECORD_LENGTH");fprintf(fheader,"%i\n\n",record_length);
0810 
0811   //(MACG) -Wshadow int byte_order = check_byte_order();
0812   byte_order = check_byte_order();
0813   write_blockname("BYTE_ORDER");fprintf(fheader,"%i\n\n",byte_order);
0814 
0815   write_blockname("ORIG_HISTORIES");
0816   if( orig_histories == 0) printf(
0817 
0818      "\n The number of primary particles (ORIG_HISTORIES) is zero in the HEADER !\n");
0819 
0820 
0821   fprintf(fheader,"%llu\n\n",orig_histories);
0822 
0823   write_blockname("PARTICLES");fprintf(fheader,"%llu\n\n",nParticles);
0824 
0825   if(particle_number[0]>0) {
0826         write_blockname("PHOTONS");fprintf(fheader,"%llu\n\n",particle_number[0]);}
0827   if(particle_number[1]>0) {
0828         write_blockname("ELECTRONS");fprintf(fheader,"%llu\n\n",particle_number[1]);}
0829   if(particle_number[2]>0) {
0830         write_blockname("POSITRONS");fprintf(fheader,"%llu\n\n",particle_number[2]);}
0831   if(particle_number[3]>0) {
0832         write_blockname("NEUTRONS");fprintf(fheader,"%llu\n\n",particle_number[3]);}
0833   if(particle_number[4]>0) {
0834         write_blockname("PROTONS");  fprintf(fheader,"%llu\n\n",particle_number[4]);}
0835 
0836   write_blockname("TRANSPORT_PARAMETERS"); fprintf(fheader,"\n");
0837 
0838   // 3. Mandatory additional information
0839   write_blockname("MACHINE_TYPE");  fprintf(fheader,"\n");
0840 
0841   write_blockname("MONTE_CARLO_CODE_VERSION"); fprintf(fheader,"\n");
0842 
0843   write_blockname("GLOBAL_PHOTON_ENERGY_CUTOFF");
0844   fprintf(fheader," %8.5f \n",global_photon_energy_cutoff);
0845 
0846   write_blockname("GLOBAL_PARTICLE_ENERGY_CUTOFF");
0847   fprintf(fheader," %8.5f \n",global_particle_energy_cutoff);
0848 
0849   write_blockname("COORDINATE_SYSTEM_DESCRIPTION"); fprintf(fheader,"\n");
0850 
0851   // 4. Optional information
0852   fprintf(fheader,"//  OPTIONAL INFORMATION\n\n");
0853 
0854   write_blockname("BEAM_NAME"); fprintf(fheader,"\n");
0855 
0856   write_blockname("FIELD_SIZE"); fprintf(fheader,"\n");
0857 
0858   write_blockname("NOMINAL_SSD"); fprintf(fheader,"\n");
0859 
0860   write_blockname("MC_INPUT_FILENAME"); fprintf(fheader,"\n");
0861 
0862   write_blockname("VARIANCE_REDUCTION_TECHNIQUES"); fprintf(fheader,"\n");
0863 
0864   write_blockname("INITIAL_SOURCE_DESCRIPTION"); fprintf(fheader,"\n");
0865 
0866   write_blockname("PUBLISHED_REFERENCE"); fprintf(fheader,"\n");
0867 
0868   write_blockname("AUTHORS"); fprintf(fheader,"\n");
0869 
0870   write_blockname("INSTITUTION"); fprintf(fheader,"\n");
0871 
0872   write_blockname("LINK_VALIDATION"); fprintf(fheader,"\n");
0873 
0874   write_blockname("ADDITIONAL_NOTES");
0875   fprintf(fheader,"%s\n","This is IAEA header as defined in the technical ");
0876   fprintf(fheader,"%s\n","report IAEA(NDS)-0484, Vienna, 2006");
0877 
0878   fprintf(fheader,"\n");
0879   // 5. Statistical information
0880   write_blockname("STATISTICAL_INFORMATION_PARTICLES");
0881   fprintf(fheader,
0882 
0883   "//        Weight        Wmin       Wmax       <E>         Emin         Emax    Particle\n");
0884   double eaver; char buffer[15];
0885   for(i=0;i<MAX_NUM_PARTICLES;i++)
0886   {
0887         if( particle_number[i] == 0 ) continue;
0888 
0889         switch (i)
0890         {
0891               case 0:
0892                     strcpy(buffer," PHOTONS");
0893                     break;
0894               case 1:
0895                     strcpy(buffer," ELECTRONS");
0896                     break;
0897               case 2:
0898                     strcpy(buffer," POSITRONS");
0899                     break;
0900               case 3:
0901                     strcpy(buffer," NEUTRONS");
0902                     break;
0903               case 4:
0904                     strcpy(buffer," PROTONS");
0905                     break;
0906         }
0907 
0908         eaver = 0.;
0909         if(sumParticleWeight[i]>0)
0910             eaver = averageKineticEnergy[i]/sumParticleWeight[i];
0911         fprintf(fheader,"  %15.6G %10.4G %10.4G %10.4G    %10.4G  %10.4G  %s\n",
0912         sumParticleWeight[i],minimumWeight[i],maximumWeight[i],
0913         eaver,minimumKineticEnergy[i],maximumKineticEnergy[i],buffer);
0914   }
0915   fprintf(fheader,"\n");
0916 
0917   write_blockname("STATISTICAL_INFORMATION_GEOMETRY");
0918   if(record_contents[0] == 1) fprintf(fheader," %G  %G\n",minimumX,maximumX);
0919   if(record_contents[1] == 1) fprintf(fheader," %G  %G\n",minimumY,maximumY);
0920   if(record_contents[2] == 1) fprintf(fheader," %G  %G\n\n",minimumZ,maximumZ);
0921 
0922   return(OK);
0923 
0924 }
0925 
0926 int iaea_header_type::print_header ()
0927 {
0928 
0929     if(checksum == 0) printf("\n NEW PHASE SPACE FILE WILL BE CREATED\n");
0930 
0931     printf("\n\nIAEA_INDEX: %i\n",iaea_index);
0932     printf("TITLE: %s \n",title);
0933 
0934   // ******************************************************************************
0935   // 1. PHSP format
0936 
0937       /*********************************************/
0938     if(file_type == 0) printf("FILE TYPE: PHASESPACE \n");
0939     if(file_type == 1) printf("FILE TYPE: GENERATOR \n");
0940 
0941     if(checksum>0) printf("CHECKSUM: %llu\n",checksum);
0942 
0943     printf("RECORD LENGTH: %i\n",record_length);
0944 
0945       if(byte_order > 0) printf("BYTE ORDER: %i\n",byte_order);
0946 
0947     int i;
0948     printf("\nRECORD_CONTENTS:\n");
0949     for (i=0;i<7;i++)
0950     {
0951             if(record_contents[i] == 0)
0952             printf(" // Variable %1i is constant\n",i+1);
0953     }
0954 
0955     if(record_contents[7] > 0)
0956         printf(" // %1i extra FLOAT variable(s) defined\n",record_contents[7]);
0957 
0958     for(i=0;i<record_contents[7];i++)
0959     {
0960       if(extrafloat_contents[i] == 0) printf(
0961         " // Generic float variable stored in the extrafloat array [%2i] \n",i);
0962       if(extrafloat_contents[i] == 1) printf(
0963         " // XLAST variable stored in the extrafloat array [%2i] \n",i);
0964       if(extrafloat_contents[i] == 2) printf(
0965         " // YLAST variable stored in the extrafloat array [%2i] \n",i);
0966       if(extrafloat_contents[i] == 3) printf(
0967         " // ZLAST variable stored in the extrafloat array [%2i] \n",i);
0968     }
0969 
0970     if(record_contents[8] > 0)
0971         printf(" // %1i extra LONG variable(s) defined\n",record_contents[8]);
0972 
0973     for(i=0;i<record_contents[8];i++)
0974     {
0975       if(extralong_contents[i] == 0) printf(
0976         " // Generic integer variable stored in the extralong array [%2i] \n",i);
0977       if(extralong_contents[i] == 1) printf(
0978         " // Incremental history number stored in the extralong array [%2i] \n",i);
0979       if(extralong_contents[i] == 2) printf(
0980         " // LATCH EGS variable stored in the extralong array [%2i] \n",i);
0981       if(extralong_contents[i] == 3) printf(
0982         " // ILB5 PENELOPE variable stored in the extralong array [%2i] \n",i);
0983       if(extralong_contents[i] == 4) printf(
0984         " // ILB4 PENELOPE variable stored in the extralong array [%2i] \n",i);
0985       if(extralong_contents[i] == 5) printf(
0986         " // ILB3 PENELOPE variable stored in the extralong array [%2i] \n",i);
0987       if(extralong_contents[i] == 6) printf(
0988         " // ILB2 PENELOPE variable stored in the extralong array [%2i] \n",i);
0989       if(extralong_contents[i] == 7) printf(
0990         " // ILB1 PENELOPE variable stored in the extralong array [%2i] \n",i);
0991     }
0992 
0993     /*********************************************/
0994     printf("\nRECORD_CONSTANT:\n");
0995     for (i=0;i<7;i++)
0996     {
0997         if(record_contents[i] > 0) continue;
0998         printf(" %8.4f // Constant variable # %1i\n",record_constant[i],i+1);
0999     };
1000     printf("\n");
1001 // ******************************************************************************
1002 // 2. Mandatory description of the phsp
1003 
1004       if( strncmp(coordinate_system_description,"                ",15) > 0 )
1005           printf("\nCOORDINATE_SYSTEM_DESCRIPTION: \n%s\n",
1006             coordinate_system_description);
1007 
1008       if(file_type == 1)
1009       {
1010             // For event generators
1011           printf("INPUT FILE for event generator: %s \n",input_file_for_event_generator);
1012             return OK;
1013       }
1014       printf("\n");
1015       printf("NUMBER OF PRIMARY PARTICLES: %llu \n",orig_histories);
1016 
1017       printf("PARTICLES: %llu \n",nParticles);
1018 
1019       if(particle_number[0] > 0) printf("PHOTONS: %llu \n",particle_number[0]);
1020       if(particle_number[1] > 0) printf("ELECTRONS: %llu \n",particle_number[1]);
1021       if(particle_number[2] > 0) printf("POSITRONS: %llu \n",particle_number[2]);
1022       if(particle_number[3] > 0) printf("NEUTRONS: %llu \n",particle_number[3]);
1023       if(particle_number[4] > 0) printf("PROTONS: %llu \n",particle_number[4]);
1024       printf("\n");
1025 
1026 // ******************************************************************************
1027 // 3. Mandatory additional information
1028       /*********************************************/
1029       if( strncmp(machine_type,"                ",15) > 0 )
1030         printf("MACHINE_TYPE: %s\n",machine_type);
1031 
1032       if( strncmp(MC_code_and_version,"                ",15) > 0 )
1033         printf("MONTE_CARLO_CODE_VERSION: %s \n",MC_code_and_version);
1034 
1035       printf("GLOBAL_PHOTON_ENERGY_CUTOFF: %8.5f \n",global_photon_energy_cutoff);
1036       printf("GLOBAL_PARTICLE_ENERGY_CUTOFF: %8.5f \n",global_particle_energy_cutoff);
1037       printf("\n");
1038 
1039       if( strncmp(transport_parameters,"                ",15) > 0 )
1040         printf("\nTRANSPORT_PARAMETERS:\n%s\n",transport_parameters);
1041 
1042 // ******************************************************************************
1043 // 4. Optional description
1044       if( strncmp(beam_name,"                ",15) > 0 )
1045         printf("BEAM_NAME: %s\n",beam_name);
1046       if( strncmp(field_size,"                ",15) > 0 )
1047         printf("FIELD_SIZE: %s\n",field_size);
1048       if( strncmp(nominal_SSD,"                ",15) > 0 )
1049         printf("NOMINAL_SSD: %s\n",nominal_SSD);
1050       if( strncmp(variance_reduction_techniques,"                ",15) > 0 )
1051         printf("VARIANCE_REDUCTION_TECHNIQUES:\n%s\n",
1052         variance_reduction_techniques);
1053       if( strncmp(initial_source_description,"                ",15) > 0 )
1054         printf("INITIAL_SOURCE_DESCRIPTION: \n%s\n",
1055         initial_source_description);
1056 
1057       // Documentation sub-section
1058       /*********************************************/
1059       if( strncmp(MC_input_filename,"                ",15) > 0 )
1060         printf("MC_INPUT_FILENAME: %s\n",MC_input_filename);
1061       if( strncmp(published_reference,"                ",15) > 0 )
1062         printf("PUBLISHED_REFERENCE: \n%s\n",published_reference);
1063       if( strncmp(authors,"                ",15) > 0 )
1064         printf("AUTHORS: \n%s\n",authors);
1065       if( strncmp(institution,"                ",15) > 0 )
1066         printf("INSTITUTION: \n%s\n",institution);
1067       if( strncmp(link_validation,"                ",15) > 0 )
1068         printf("LINK_VALIDATION: \n%s\n",link_validation);
1069       if( strncmp(additional_notes,"                ",15) > 0 )
1070         printf("ADDITIONAL_NOTES: \n%s\n",additional_notes);
1071 
1072 // ******************************************************************************
1073 // 5. Optional statistical information
1074 
1075       print_statistics();
1076 
1077     return(OK);
1078 }