File indexing completed on 2026-04-18 07:41:08
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046
0047
0048
0049
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
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);
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
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)
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)
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);
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);
0200
0201
0202 for(int itmp=0;itmp<MAX_NUM_PARTICLES;itmp++) particle_number[itmp]=0;
0203 IAEA_I64 npart;
0204
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
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
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
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
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
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;
0314 fbuff[icnt++] = (float)atof(line+index1);
0315 if(advance(line, &index0, len) != OK) break;
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
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;
0345 fbuff[icnt++] = (float)atof(line+index1);
0346 if(advance(line, &index0, len) != OK) break;
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
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
0384 endptr = strchr(begptr,SEGMENT_END_TOKEN) ;
0385 if( endptr != NULL ) {
0386 *endptr = '\0' ;
0387
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,"");
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;
0450 for(i=0;i<8;i++) record_length += record_contents[i]*sizeof(float);
0451 record_length -= 4;
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;
0500 for(i=0;i<8;i++) record_length += record_contents[i]*sizeof(float);
0501 record_length -= 4;
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
0683 float ftest=1.0f;
0684 char *pf = (char *) &ftest;
0685
0686 if(pf[0] == 0 && pf[3] != 0)
0687 {
0688
0689 return(LITTLE_ENDIAN);
0690 }else if(pf[0] != 0 && pf[3] == 0)
0691 {
0692
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");
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
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
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
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
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
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
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
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
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
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
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
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
1074
1075 print_statistics();
1076
1077 return(OK);
1078 }