1/******************************************************************************* 2* 3* McStas, neutron ray-tracing package 4* Copyright(C) 2007 Risoe National Laboratory. 5* 6* %I 7* Written by: Mads Bertelsen 8* Date: 20.08.15 9* Version: $Revision: 0.1 $ 10* Origin: Svanevej 19 11* 12* A sample component to separate geometry and phsysics 13* 14* %D 15* Part of the Union components, a set of components that work together and thus 16* sperates geometry and physics within McStas. 17* The use of this component requires other components to be used. 18* 19* 1) One specifies a number of processes using process components 20* 2) These are gathered into material definitions using Union_make_material 21* 3) Geometries are placed using Union_box/cylinder/sphere, assigned a material 22* 4) A Union_master component placed after all of the above 23* 24* Only in step 4 will any simulation happen, and per default all geometries 25* defined before this master, but after the previous will be simulated here. 26* 27* There is a dedicated manual available for the Union_components 28* 29* This logger logs a 2D projection of the position of each scattering in the lab 30* frame. It does this in n3 space slices, so one can get a full 3D histogram of 31* scattering positions. 32* 33* A logger will log something for scattering events happening to certain volumes, 34* which are specified in the target_geometry string. By leaving it blank, all 35* geometries are logged, even the ones not defined at this point in the 36* instrument file. If a list og target_geometries is selected, one can further 37* narrow the events logged by providing a list of process names in target_process 38* which need to correspond with names of defined Union_process components. 39* 40* To use the logger_conditional_extend function, set it to some integer value n 41* and make and extend section to the master component that runs the geometry. 42* In this extend function, logger_conditional_extend[n] is 1 if the conditional 43* stack evaluated to true, 0 if not. This way one can check what rays is logged 44* using regular McStas monitors. Only works if a conditional is applied to this 45* logger. 46* 47* 48* %P 49* INPUT PARAMETERS: 50* D_direction_1: [string] Direction for first axis ("x", "y" or "z") 51* D1_max: [1] histogram boundery, max position value for first axis 52* D1_min: [1] histogram boundery, min position value for first axis 53* n1: [1] number of bins for first axis 54* D_direction_2: [string] Direction for second axis ("x", "y" or "z") 55* D2_max: [1] histogram boundery, max position value for second axis 56* D2_min: [1] histogram boundery, min position value for second axis 57* n2: [1] number of bins for second axis 58* D_direction_3: [string] Direction for second axis ("x", "y" or "z") 59* D3_max: [1] histogram boundery, max position value for third axis 60* D3_min: [1] histogram boundery, min position value for third axis 61* n3: [1] number of bins for third axis 62* target_geometry: [string] Comma seperated list of geometry names that will be logged, leave empty for all volumes (even not defined yet) 63* target_process: [string] Comma seperated names of physical processes, if volumes are selected, one can select Union_process names 64* order_total: [1] Only log rays that scatter for the n'th time, 0 for all orders 65* order_volume: [1] Only log rays that scatter for the n'th time in the same geometry 66* order_volume_process [1] Only log rays that scatter for the n'th time in the same geometry, using the same process 67* logger_conditional_extend_index [1] If a conditional is used with this logger, the result of each conditional calculation can be made available in extend as a array called "logger_conditional_extend", and one would then acces logger_conditional_extend[n] if logger_conditional_extend_index is set to n 68* 69* OUTPUT PARAMETERS: 70* 71* GLOBAL PARAMETERS: 72* 73* %L 74* 75* %E 76******************************************************************************/ 77 78DEFINE COMPONENT Union_logger_3D_space 79DEFINITION PARAMETERS (n1=90, n2=90, n3=10) 80SETTING PARAMETERS(string target_geometry="NULL",string target_process="NULL",D1_min=-1,D1_max=1,D2_min=-1,D2_max=1,D3_min=-1,D3_max=1,string D_direction_1="x", string D_direction_2="z", string D_direction_3="z",string filename="NULL", order_total=0,order_volume=0,order_volume_process=0,logger_conditional_extend_index=-1) 81OUTPUT PARAMETERS (loop_index,accepted_processes,accepted_volumes,logger_list_element,this_logger,this_storage,loggers_on_target_volume,target_volume,logger_conditional_extend_index,part_filename,number_string,storage_N,storage_p,storage_2p) 82 83/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ 84 85SHARE 86%{ 87#ifndef Union 88#define Union $Revision: 0.8 $ 89 90#include "Union_functions.c" 91#include "Union_initialization.c" 92 93#endif 94 95struct temp_3DS_data_element_struct { 96 int index_1; 97 int index_2; 98 int index_3; 99 double weight; 100}; 101 102struct temp_3DS_data_struct { 103 int num_elements; 104 int allocated_elements; 105 struct temp_3DS_data_element_struct *elements; 106}; 107 108struct a_3DS_storage_struct { 109 struct Detector_3D_struct Detector_3D; 110 struct temp_3DS_data_struct temp_3DS_data; 111 int dim_1_choice; 112 int dim_2_choice; 113 int dim_3_choice; 114 int order; 115 int order_in_this_volume; 116 int order_process_in_this_volume; 117 118 Coords position; 119 Rotation rotation; 120 Rotation t_rotation; 121}; 122 123// record_to_temp 124// Would be nice if x y z, k_new and k_old were all coords 125void record_to_temp_3DS(Coords *position, double *k_new, double *k_old, double p, double p_old, double time, int scattered_in_this_volume, int scattered_in_this_volume_by_this_process, int total_number_of_scattering_events, struct logger_struct *logger, struct logger_with_data_struct *logger_with_data_array) { 126 127 struct a_3DS_storage_struct *storage; 128 storage = logger->data_union.p_3DS_storage; 129 130 int add_point = 1; 131 132 if (storage->order != 0) { 133 if (storage->order - 1 == total_number_of_scattering_events) 134 add_point = 1; 135 else 136 add_point = 0; 137 } 138 139 if (storage->order_in_this_volume != 0) { 140 if (storage->order_in_this_volume - 1 == scattered_in_this_volume) 141 add_point = 1; 142 else 143 add_point = 0; 144 } 145 146 if (storage->order_process_in_this_volume != 0) { 147 if (storage->order_process_in_this_volume - 1 == scattered_in_this_volume_by_this_process) 148 add_point = 1; 149 else 150 add_point = 0; 151 } 152 153 if (add_point == 1) { 154 155 double p1,p2,p3; 156 157 // dim_1_choice = 0 for x, 1 for y, 2 for z. Done in initialize from input. "x" "y" "z". 158 if (storage->dim_1_choice == 0) 159 p1 = position->x - storage->position.x; 160 else if (storage->dim_1_choice == 1) 161 p1 = position->y - storage->position.y; 162 else if (storage->dim_1_choice == 2) 163 p1 = position->z - storage->position.z; 164 165 if (storage->dim_2_choice == 0) 166 p2 = position->x - storage->position.x; 167 else if (storage->dim_2_choice == 1) 168 p2 = position->y - storage->position.y; 169 else if (storage->dim_2_choice == 2) 170 p2 = position->z - storage->position.z; 171 172 if (storage->dim_3_choice == 0) 173 p2 = position->x - storage->position.x; 174 else if (storage->dim_3_choice == 1) 175 p2 = position->y - storage->position.y; 176 else if (storage->dim_3_choice == 2) 177 p2 = position->z - storage->position.z; 178 179 int i,j,k; 180 181 // Find bin in histogram 182 if (p1>storage->Detector_3D.D1min && p1<storage->Detector_3D.D1max && p2>storage->Detector_3D.D2min && p2<storage->Detector_3D.D2max && p3>storage->Detector_3D.D3min && p3<storage->Detector_3D.D3max) { 183 i = floor((p1 - storage->Detector_3D.D1min)*storage->Detector_3D.bins_1/(storage->Detector_3D.D1max - storage->Detector_3D.D1min)); 184 j = floor((p2 - storage->Detector_3D.D2min)*storage->Detector_3D.bins_2/(storage->Detector_3D.D2max - storage->Detector_3D.D2min)); 185 k = floor((p3 - storage->Detector_3D.D3min)*storage->Detector_3D.bins_3/(storage->Detector_3D.D3max - storage->Detector_3D.D3min)); 186 187 188 // Save bin in histogram to temp (may need to allocate more memory) 189 int index; 190 //printf("number of data points used: %d space allocated for %d data points. \n",storage->temp_3DS_data.num_elements,storage->temp_3DS_data.allocated_elements); 191 192 if (storage->temp_3DS_data.num_elements < storage->temp_3DS_data.allocated_elements) { 193 storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_1 = k; 194 storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_2 = i; 195 storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_3 = j; 196 storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements++].weight = p; 197 } else { 198 // No more space, need to allocate a larger buffer for this logger. Wish I had generics. 199 200 // copy current data to temp 201 struct temp_3DS_data_struct temporary_storage; 202 temporary_storage.num_elements = storage->temp_3DS_data.num_elements; 203 temporary_storage.elements = malloc(temporary_storage.num_elements*sizeof(struct temp_3DS_data_element_struct)); 204 205 for (index=0;index<storage->temp_3DS_data.num_elements;index++) { 206 temporary_storage.elements[index].index_1 = storage->temp_3DS_data.elements[index].index_1; 207 temporary_storage.elements[index].index_2 = storage->temp_3DS_data.elements[index].index_2; 208 temporary_storage.elements[index].index_3 = storage->temp_3DS_data.elements[index].index_3; 209 temporary_storage.elements[index].weight = storage->temp_3DS_data.elements[index].weight; 210 } 211 212 // free current data 213 free(storage->temp_3DS_data.elements); 214 215 // allocate larger array (10 larger) 216 storage->temp_3DS_data.allocated_elements = 10 + storage->temp_3DS_data.num_elements; 217 storage->temp_3DS_data.elements = malloc(storage->temp_3DS_data.allocated_elements*sizeof(struct temp_3DS_data_element_struct)); 218 219 // copy back from temp 220 for (index=0;index<storage->temp_3DS_data.num_elements;index++) { 221 storage->temp_3DS_data.elements[index].index_1 = temporary_storage.elements[index].index_1; 222 storage->temp_3DS_data.elements[index].index_2 = temporary_storage.elements[index].index_2; 223 storage->temp_3DS_data.elements[index].index_3 = temporary_storage.elements[index].index_3; 224 storage->temp_3DS_data.elements[index].weight = temporary_storage.elements[index].weight; 225 } 226 227 // free temporary data 228 free(temporary_storage.elements); 229 230 // add new data point 231 storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_1 = k; 232 storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_2 = i; 233 storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_3 = j; 234 storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements++].weight = p; 235 } 236 237 // If this is the first time this ray is being recorded in this logger, add it to the list of loggers that write to temp and may get it moved to perm 238 if (storage->temp_3DS_data.num_elements == 1) 239 add_to_logger_with_data(logger_with_data_array,logger); 240 } 241 } 242 243} 244 245// clear_temp 246void clear_temp_3DS(union logger_data_union *data_union) { 247 data_union->p_3DS_storage->temp_3DS_data.num_elements = 0; 248} 249 250// record_to_perm 251void record_to_perm_3DS(Coords *position, double *k_new, double *k_old, double p, double p_old, double time, int scattered_in_this_volume, int scattered_in_this_volume_by_this_process, int total_number_of_scattering_events, struct logger_struct *logger, struct logger_with_data_struct *logger_with_data_array) { 252 253 //printf("In record to permanent \n"); 254 struct a_3DS_storage_struct *storage; 255 storage = logger->data_union.p_3DS_storage; 256 257 int add_point = 1; 258 259 if (storage->order != 0) { 260 if (storage->order - 1 == total_number_of_scattering_events) 261 add_point = 1; 262 else 263 add_point = 0; 264 } 265 266 if (storage->order_in_this_volume != 0) { 267 if (storage->order_in_this_volume - 1 == scattered_in_this_volume) 268 add_point = 1; 269 else 270 add_point = 0; 271 } 272 273 if (storage->order_process_in_this_volume != 0) { 274 if (storage->order_process_in_this_volume - 1 == scattered_in_this_volume_by_this_process) 275 add_point = 1; 276 else 277 add_point = 0; 278 } 279 280 if (add_point == 1) { 281 //printf("storage was set \n"); 282 double p1,p2,p3; 283 284 // dim_1_choice = 0 for x, 1 for y, 2 for z. Done in initialize from input. "x" "y" "z". 285 if (storage->dim_1_choice == 0) 286 p1 = position->x - storage->position.x; 287 else if (storage->dim_1_choice == 1) 288 p1 = position->y - storage->position.y; 289 else if (storage->dim_1_choice == 2) 290 p1 = position->z - storage->position.z; 291 292 if (storage->dim_2_choice == 0) 293 p2 = position->x - storage->position.x; 294 else if (storage->dim_2_choice == 1) 295 p2 = position->y - storage->position.y; 296 else if (storage->dim_2_choice == 2) 297 p2 = position->z - storage->position.z; 298 299 if (storage->dim_3_choice == 0) 300 p3 = position->x - storage->position.x; 301 else if (storage->dim_3_choice == 1) 302 p3 = position->y - storage->position.y; 303 else if (storage->dim_3_choice == 2) 304 p3 = position->z - storage->position.z; 305 306 int i,j,k; 307 308 // Find bin in histogram 309 if (p1>storage->Detector_3D.D1min && p1<storage->Detector_3D.D1max && p2>storage->Detector_3D.D2min && p2<storage->Detector_3D.D2max && p3>storage->Detector_3D.D3min && p3<storage->Detector_3D.D3max) { 310 i = floor((p1 - storage->Detector_3D.D1min)*storage->Detector_3D.bins_1/(storage->Detector_3D.D1max - storage->Detector_3D.D1min)); 311 j = floor((p2 - storage->Detector_3D.D2min)*storage->Detector_3D.bins_2/(storage->Detector_3D.D2max - storage->Detector_3D.D2min)); 312 k = floor((p3 - storage->Detector_3D.D3min)*storage->Detector_3D.bins_3/(storage->Detector_3D.D3max - storage->Detector_3D.D3min)); 313 314 /* 315 printf("(p1,p2,p3) = (%f,%f,%f)\n",p1,p2,p3); 316 printf("limits1 = [%f,%f] \n",storage->Detector_3D.D1min,storage->Detector_3D.D1max); 317 printf("limits2 = [%f,%f] \n",storage->Detector_3D.D2min,storage->Detector_3D.D2max); 318 printf("limits3 = [%f,%f] \n",storage->Detector_3D.D3min,storage->Detector_3D.D3max); 319 printf("n bins = [%d,%d,%d] \n",round(storage->Detector_3D.bins_1),round(storage->Detector_3D.bins_2),round(storage->Detector_3D.bins_3)); 320 printf("Added to statistics for monitor [%d] [%d] [%d] \n",i,j,k); 321 printf("indicies found\n"); 322 */ 323 324 /* 325 storage->Detector_3D.Array_N[i][j][k]++; 326 storage->Detector_3D.Array_p[i][j][k] += p; 327 storage->Detector_3D.Array_p2[i][j][k] += p*p; 328 */ 329 330 331 // because of the order in which the elements are laid out in memory, the k index must be first. 332 storage->Detector_3D.Array_N[k][i][j]++; 333 storage->Detector_3D.Array_p[k][i][j] += p; 334 storage->Detector_3D.Array_p2[k][i][j] += p*p; 335 336 337 } 338 } 339} 340 341// write_temp_to_perm 342void write_temp_to_perm_3DS(union logger_data_union *data_union) { 343 344 struct a_3DS_storage_struct *storage; 345 storage = data_union->p_3DS_storage; 346 347 int index; 348 // Add all data points to the historgram, they are saved as index / weight combinations 349 for (index=0;index<storage->temp_3DS_data.num_elements;index++) { 350 storage->Detector_3D.Array_N[storage->temp_3DS_data.elements[index].index_1][storage->temp_3DS_data.elements[index].index_2][storage->temp_3DS_data.elements[index].index_3]++; 351 352 storage->Detector_3D.Array_p[storage->temp_3DS_data.elements[index].index_1][storage->temp_3DS_data.elements[index].index_2][storage->temp_3DS_data.elements[index].index_3] += storage->temp_3DS_data.elements[index].weight; 353 354 storage->Detector_3D.Array_p2[storage->temp_3DS_data.elements[index].index_1][storage->temp_3DS_data.elements[index].index_2][storage->temp_3DS_data.elements[index].index_3] += storage->temp_3DS_data.elements[index].weight*storage->temp_3DS_data.elements[index].weight; 355 } 356 clear_temp_3DS(data_union); 357} 358 359void write_temp_to_perm_final_p_3DS(union logger_data_union *data_union, double final_weight) { 360 361 struct a_3DS_storage_struct *storage; 362 storage = data_union->p_3DS_storage; 363 364 int index; 365 // Add all data points to the historgram, they are saved as index / weight combinations 366 for (index=0;index<storage->temp_3DS_data.num_elements;index++) { 367 storage->Detector_3D.Array_N[storage->temp_3DS_data.elements[index].index_1][storage->temp_3DS_data.elements[index].index_2][storage->temp_3DS_data.elements[index].index_3]++; 368 369 storage->Detector_3D.Array_p[storage->temp_3DS_data.elements[index].index_1][storage->temp_3DS_data.elements[index].index_2][storage->temp_3DS_data.elements[index].index_3] += final_weight; 370 371 storage->Detector_3D.Array_p2[storage->temp_3DS_data.elements[index].index_1][storage->temp_3DS_data.elements[index].index_2][storage->temp_3DS_data.elements[index].index_3] += final_weight*final_weight; 372 } 373 clear_temp_3DS(data_union); 374} 375 376// Only need to define linking function for loggers once. 377#ifndef UNION_LOGGER 378#define UNION_LOGGER Dummy 379// Linking function for loggers, finds the indicies of the specified geometries on the global_geometry_list 380void manual_linking_function_logger_volumes(char *input_string, struct pointer_to_global_geometry_list *global_geometry_list, struct pointer_to_1d_int_list *accepted_volumes, char *component_name) { 381 // Need to check a input_string of text for an occurance of name. If it is in the inputstring, yes return 1, otherwise 0. 382 char *token; 383 int loop_index; 384 char local_string[512]; 385 386 strcpy(local_string,input_string); 387 // get the first token 388 token = strtok(local_string,","); 389 390 // walk through other tokens 391 while( token != NULL ) 392 { 393 //printf( " %s\n", token ); 394 for (loop_index=0;loop_index<global_geometry_list->num_elements;loop_index++) { 395 if (strcmp(token,global_geometry_list->elements[loop_index].name) == 0) { 396 add_element_to_int_list(accepted_volumes,loop_index); 397 break; 398 } 399 400 if (loop_index == global_geometry_list->num_elements - 1) { 401 // All possible geometry names have been looked through, and the break was not executed. 402 // Alert the user to this problem by showing the geometry name that was not found and the currently available geometires 403 printf("\n"); 404 printf("ERROR: The target_geometry string \"%s\" in Union logger component \"%s\" had an entry that did not match a specified geometry. \n",input_string,component_name); 405 printf(" The unrecoignized geometry name was: \"%s\" \n",token); 406 printf(" The geometries available at this point (need to be defined before the logger): \n"); 407 for (loop_index=0;loop_index<global_geometry_list->num_elements;loop_index++) 408 printf(" %s\n",global_geometry_list->elements[loop_index].name); 409 exit(1); 410 } 411 } 412 413 // Updates the token 414 token = strtok(NULL,","); 415 } 416} 417 418void manual_linking_function_logger_processes(char *input_string, struct physics_struct *p_physics, struct pointer_to_1d_int_list *accepted_processes, char *component_name, char *Volume_name) { 419 // Need to check a input_string of text for an occurance of name. If it is in the inputstring, yes return 1, otherwise 0. 420 char *token; 421 int loop_index; 422 char local_string[256]; 423 424 strcpy(local_string,input_string); 425 // get the first token 426 token = strtok(local_string,","); 427 428 // walk through other tokens 429 while( token != NULL ) 430 { 431 //printf( " %s\n", token ); 432 for (loop_index=0;loop_index<p_physics->number_of_processes;loop_index++) { 433 if (strcmp(token,p_physics->p_scattering_array[loop_index].name) == 0) { 434 add_element_to_int_list(accepted_processes,loop_index); 435 break; 436 } 437 438 if (loop_index == p_physics->number_of_processes - 1) { 439 // All possible process names have been looked through, and the break was not executed. 440 // Alert the user to this problem by showing the process name that was not found and the currently available processes 441 printf("\n"); 442 printf("ERROR: The target process string \"%s\" in Union logger \"%s\" had an entry that did not match a specified process in assosiated volume \"%s\". \n",input_string,component_name,Volume_name); 443 printf(" The unrecoignized process name was: \"%s\" \n",token); 444 printf(" The processes defined in material \"%s\" of which Volume \"%s\" is made: \n",p_physics->name,Volume_name); 445 for (loop_index=0;loop_index<p_physics->number_of_processes;loop_index++) 446 printf(" %s\n",p_physics->p_scattering_array[loop_index].name); 447 exit(1); 448 } 449 } 450 451 // Updates the token 452 token = strtok(NULL,","); 453 } 454} 455#endif 456 457double*** allocate3Ddouble_3DS(int count_x, int count_y, int count_z, double *storage) { 458 //double *storage = malloc(count_x * count_y * count_z * sizeof(double)); 459 storage = malloc(count_x * count_y * count_z * sizeof(double)); 460 double *alloc = storage; 461 double ***x; 462 int i,j; 463 x = malloc(count_x * sizeof(*x)); 464 for (i = 0; i < count_x; i++) { 465 x[i] = malloc(count_y * sizeof(**x)); 466 for (j = 0; j < count_y; j++) { 467 x[i][j] = alloc; 468 alloc += count_z; 469 } 470 } 471 return x; 472 } 473 474void free3Ddouble_3DS(double*** ptr_array, int count_x, double *storage) { 475 if (!ptr_array) return; 476 int x; 477 478 free(storage); 479 for (x=0;x<count_x;x++) { 480 free(ptr_array[x]); 481 } 482 free(ptr_array); 483} 484 485%} 486 487DECLARE 488%{ 489// From make material 490// Needed for transport to the main component 491//struct global_material_element_struct global_material_element; 492//struct physics_struct this_material; 493 494int loop_index,l1,l2,l3; 495int process_index; 496int found_process; 497int specified_processes; 498char local_string[256]; 499 500char number_string[16]; 501char part_filename[256]; 502 503// Reused for logger 504struct pointer_to_1d_int_list accepted_processes = {0,NULL}; 505 506struct global_logger_element_struct logger_list_element; 507 508struct pointer_to_1d_int_list accepted_volumes = {0,NULL}; 509 510struct logger_struct this_logger; 511struct a_3DS_storage_struct this_storage; 512 513struct loggers_struct *loggers_on_target_volume; 514struct Volume_struct *target_volume; 515 516char temp_string[2]; 517 518double *storage_N,*storage_p,*storage_2p; 519 520%} 521 522INITIALIZE 523%{ 524 525 526 // Initialize storage from input 527 if (D1_min >= D1_max) { 528 printf("ERROR, Union logger \"%s\" had D1_min >= D1_max.\n",NAME_CURRENT_COMP); 529 exit(1); 530 } 531 this_storage.Detector_3D.D1min = D1_min; 532 this_storage.Detector_3D.D1max = D1_max; 533 534 if (D2_min >= D2_max) { 535 printf("ERROR, Union logger \"%s\" had D2_min >= D2_max.\n",NAME_CURRENT_COMP); 536 exit(1); 537 } 538 this_storage.Detector_3D.D2min = D2_min; 539 this_storage.Detector_3D.D2max = D2_max; 540 541 if (D3_min >= D3_max) { 542 printf("ERROR, Union logger \"%s\" had D3_min >= D3_max.\n",NAME_CURRENT_COMP); 543 exit(1); 544 } 545 this_storage.Detector_3D.D3min = D3_min; 546 this_storage.Detector_3D.D3max = D3_max; 547 548 if (n1 <= 0) { 549 printf("ERROR, Union logger \"%s\" had n1 <= 0.\n",NAME_CURRENT_COMP); 550 exit(1); 551 } 552 this_storage.Detector_3D.bins_1 = n1; 553 554 if (n2 <= 0) { 555 printf("ERROR, Union logger \"%s\" had n2 <= 0.\n",NAME_CURRENT_COMP); 556 exit(1); 557 } 558 this_storage.Detector_3D.bins_2 = n2; 559 560 if (n3 <= 0) { 561 printf("ERROR, Union logger \"%s\" had n3 <= 0.\n",NAME_CURRENT_COMP); 562 exit(1); 563 } 564 this_storage.Detector_3D.bins_3 = n3; 565 566 //printf("past input sanitation \n"); 567 568 569 printf("Allocating 3D arrays \n"); 570 // Remember to take special care when deallocating this array, use free3Ddouble 571 572 /* 573 this_storage.Detector_3D.Array_N = allocate3Ddouble_3DS(n1,n2,n3); // Here the n1 double is cast to an int 574 this_storage.Detector_3D.Array_p = allocate3Ddouble_3DS(n1,n2,n3); 575 this_storage.Detector_3D.Array_p2 = allocate3Ddouble_3DS(n1,n2,n3); 576 */ 577 // D3 is in poisition 1, because that gives continous memory in XY plane for easy plotting 578 579 this_storage.Detector_3D.Array_N = allocate3Ddouble_3DS(n3,n1,n2,storage_N); // Here the n1 double is cast to an int 580 this_storage.Detector_3D.Array_p = allocate3Ddouble_3DS(n3,n1,n2,storage_p); 581 this_storage.Detector_3D.Array_p2 = allocate3Ddouble_3DS(n3,n1,n2,storage_2p); 582 583 /* 584 // Error in the order? 585 this_storage.Detector_3D.Array_N = allocate3Ddouble_3DS(n1,n3,n2); // Here the n1 double is cast to an int 586 this_storage.Detector_3D.Array_p = allocate3Ddouble_3DS(n1,n3,n2); 587 this_storage.Detector_3D.Array_p2 = allocate3Ddouble_3DS(n1,n3,n2); 588 */ 589 590 591 //printf("Allocated 3D arrays \n"); 592 for (l1=0;l1<n1;l1++) { //n1 is technically a double, but this works fine 593 for (l2=0;l2<n2;l2++) { 594 for (l3=0;l3<n3;l3++) { 595 this_storage.Detector_3D.Array_N[l3][l1][l2] = 0; 596 this_storage.Detector_3D.Array_p[l3][l1][l2] = 0; 597 this_storage.Detector_3D.Array_p2[l3][l1][l2] = 0; 598 599 } 600 } 601 } 602 603 //printf("Initialized 3D arrays \n"); 604 605 //printf("past 3D pointer assignment \n"); 606 607 // Input sanitation for filename apparently done in 3D_detector_out 608 609 if (strcmp(D_direction_1,"x") == 0 || strcmp(D_direction_1,"X") == 0) { 610 this_storage.dim_1_choice = 0; 611 sprintf(this_storage.Detector_3D.string_axis_1,"x [m]"); 612 sprintf(this_storage.Detector_3D.title_string,"3D position Union logger in plane x"); 613 } else if (strcmp(D_direction_1,"y") == 0 || strcmp(D_direction_1,"Y") == 0) { 614 this_storage.dim_1_choice = 1; 615 sprintf(this_storage.Detector_3D.string_axis_1,"y [m]"); 616 sprintf(this_storage.Detector_3D.title_string,"3D position Union logger in plane y"); 617 } else if (strcmp(D_direction_1,"z") == 0 || strcmp(D_direction_1,"Z") == 0) { 618 this_storage.dim_1_choice = 2; 619 sprintf(this_storage.Detector_3D.string_axis_1,"z [m]"); 620 sprintf(this_storage.Detector_3D.title_string,"3D position Union logger in plane z"); 621 } else { 622 printf("ERROR, Union logger 3DS named \"%s\" has an invalid D_direction_1 string, it should be \"x\",\"y\" or \"z\".\n",NAME_CURRENT_COMP); 623 exit(1); 624 } 625 626 if (strcmp(D_direction_2,"x") == 0 || strcmp(D_direction_2,"X") == 0) { 627 this_storage.dim_2_choice = 0; 628 sprintf(this_storage.Detector_3D.string_axis_2,"x [m]"); 629 sprintf(temp_string,"x"); 630 } else if (strcmp(D_direction_2,"y") == 0 || strcmp(D_direction_2,"Y") == 0) { 631 this_storage.dim_2_choice = 1; 632 sprintf(this_storage.Detector_3D.string_axis_2,"y [m]"); 633 sprintf(temp_string,"y"); 634 } else if (strcmp(D_direction_2,"z") == 0 || strcmp(D_direction_2,"Z") == 0) { 635 this_storage.dim_2_choice = 2; 636 sprintf(this_storage.Detector_3D.string_axis_2,"z [m]"); 637 sprintf(temp_string,"z"); 638 } else { 639 printf("ERROR, Union logger 3DS named \"%s\" has an invalid D_direction_2 string, it should be \"x\",\"y\" or \"z\".\n",NAME_CURRENT_COMP); 640 exit(1); 641 } 642 643 strcat(this_storage.Detector_3D.title_string,temp_string); // Connects the title string started in D_direction_1 part with the ending in D_direction_2 part 644 645 if (strcmp(D_direction_3,"x") == 0 || strcmp(D_direction_3,"X") == 0) { 646 this_storage.dim_3_choice = 0; 647 sprintf(this_storage.Detector_3D.string_axis_3,"x [m]"); 648 sprintf(temp_string,"x"); 649 } else if (strcmp(D_direction_3,"y") == 0 || strcmp(D_direction_3,"Y") == 0) { 650 this_storage.dim_3_choice = 1; 651 sprintf(this_storage.Detector_3D.string_axis_3,"y [m]"); 652 sprintf(temp_string,"y"); 653 } else if (strcmp(D_direction_3,"z") == 0 || strcmp(D_direction_3,"Z") == 0) { 654 this_storage.dim_3_choice = 2; 655 sprintf(this_storage.Detector_3D.string_axis_3,"z [m]"); 656 sprintf(temp_string,"z"); 657 } else { 658 printf("ERROR, Union logger 3DS named \"%s\" has an invalid D_direction_3 string, it should be \"x\",\"y\" or \"z\".\n",NAME_CURRENT_COMP); 659 exit(1); 660 } 661 662 strcat(this_storage.Detector_3D.title_string,temp_string); 663 664 sprintf(this_storage.Detector_3D.Filename,filename); 665 666 667 this_storage.order = order_total; 668 this_storage.order_in_this_volume = order_volume; 669 this_storage.order_process_in_this_volume = order_volume_process; 670 this_storage.temp_3DS_data.num_elements=0; 671 this_storage.temp_3DS_data.allocated_elements = 10; 672 this_storage.temp_3DS_data.elements = malloc(this_storage.temp_3DS_data.allocated_elements*sizeof(struct temp_3DS_data_element_struct)); 673 674 675 this_storage.position = POS_A_CURRENT_COMP; 676 add_position_pointer_to_list(&global_positions_to_transform_list,&this_storage.position); 677 678 rot_copy(this_storage.rotation,ROT_A_CURRENT_COMP); 679 add_rotation_pointer_to_list(&global_rotations_to_transform_list,&this_storage.rotation); 680 681 rot_transpose(ROT_A_CURRENT_COMP,this_storage.t_rotation); 682 add_rotation_pointer_to_list(&global_rotations_to_transform_list,&this_storage.t_rotation); 683 684 685 //printf("past direction choices sanitation \n"); 686 687 // Book keeping 688 this_logger.logger_extend_index = logger_conditional_extend_index; 689 this_logger.function_pointers.active_record_function = &record_to_perm_3DS; // Assume no conditional 690 this_logger.function_pointers.inactive_record_function = &record_to_temp_3DS; // If an assume is present, these two pointers are switched 691 692 // Temp to perm functions, and standard identifier 693 this_logger.function_pointers.select_t_to_p = 1; // 1: temp_to_perm, 2: temp_to_perm_final_p 694 this_logger.function_pointers.temp_to_perm = &write_temp_to_perm_3DS; 695 this_logger.function_pointers.temp_to_perm_final_p = &write_temp_to_perm_final_p_3DS; 696 this_logger.function_pointers.clear_temp = &clear_temp_3DS; 697 // Initializing for conditional 698 this_logger.conditional_list.num_elements = 0; 699 700 //this_logger.function_pointers.perm_to_disk = &write_perm_to_disk_3DS; //Obsolete 701 702 //printf("past this_logger function assignment \n"); 703 704 this_logger.data_union.p_3DS_storage = &this_storage; 705 706 sprintf(this_logger.name,NAME_CURRENT_COMP); 707 708 //printf("past this_logger assignment \n"); 709 710 sprintf(logger_list_element.name,NAME_CURRENT_COMP); 711 logger_list_element.component_index = INDEX_CURRENT_COMP; 712 logger_list_element.logger = &this_logger; 713 714 //printf("past logger_list_element assignment \n"); 715 716 // In order to run the logger at the right times, pointers to this logger is stored in each volume it logs, 717 // and additionally for each avaiable process. If a process is not logged, the pointer is simply not stored. 718 719 // Need to find the volumes for which the processes should have a reference to this logger 720 if (target_geometry && strlen(target_geometry) && strcmp(target_geometry,"NULL") && strcmp(target_geometry, "0")) { 721 // Certain volumes were selected, find the indicies in the global_geometry_list 722 manual_linking_function_logger_volumes(target_geometry, &global_geometry_list, &accepted_volumes, NAME_CURRENT_COMP); 723 // Add this logger to the global_specific_volumes_logger_list (so that conditionals can affect it) 724 add_element_to_logger_list(&global_specific_volumes_logger_list,logger_list_element); 725 726 for (loop_index=0;loop_index<accepted_volumes.num_elements;loop_index++) { 727 target_volume = global_geometry_list.elements[accepted_volumes.elements[loop_index]].Volume; 728 // Add an element to its logger list 729 add_initialized_logger_in_volume(&target_volume->loggers,target_volume->p_physics->number_of_processes); 730 731 if (target_process && strlen(target_process) && strcmp(target_process,"NULL") && strcmp(target_process, "0")) { 732 // Unused process pointers should point to NULL 733 for (process_index=0;process_index<target_volume->p_physics->number_of_processes;process_index++) { 734 target_volume->loggers.p_logger_volume[target_volume->loggers.num_elements-1].p_logger_process[process_index]=NULL; 735 } 736 // A target_process was set, find it within the volume structure (can be many processes) 737 manual_linking_function_logger_processes(target_process, target_volume->p_physics, &accepted_processes, NAME_CURRENT_COMP,target_volume->name); 738 for (process_index=0;process_index<accepted_processes.num_elements;process_index++) { 739 // Add pointer to this logger for all the accepted processes in this newly added loggers element 740 target_volume->loggers.p_logger_volume[target_volume->loggers.num_elements-1].p_logger_process[accepted_processes.elements[process_index]]=&this_logger; 741 } 742 } else { 743 // No target_process was set, attatch the logger to all processes 744 for (process_index=0;process_index<target_volume->p_physics->number_of_processes;process_index++) { 745 target_volume->loggers.p_logger_volume[target_volume->loggers.num_elements-1].p_logger_process[process_index]=&this_logger; 746 } 747 } 748 } 749 } else { 750 // Send to global_all_volumes_logger_list 751 // Here there is no system for selecting processes as well 752 add_element_to_logger_list(&global_all_volume_logger_list,logger_list_element); 753 } 754 755 756 %} 757 758TRACE 759%{ 760%} 761 762SAVE 763%{ 764// Write to disk 765 766for (loop_index=0;loop_index<this_storage.Detector_3D.bins_3;loop_index++) { 767 sprintf(number_string,"%d",loop_index); 768 sprintf(part_filename,"%s_%s",this_storage.Detector_3D.Filename,number_string); 769 770 DETECTOR_OUT_2D( 771 this_storage.Detector_3D.title_string, 772 this_storage.Detector_3D.string_axis_1, 773 this_storage.Detector_3D.string_axis_2, 774 this_storage.Detector_3D.D1min, this_storage.Detector_3D.D1max, 775 this_storage.Detector_3D.D2min, this_storage.Detector_3D.D2max, 776 this_storage.Detector_3D.bins_1, this_storage.Detector_3D.bins_2, 777 &this_storage.Detector_3D.Array_N[loop_index][0][0], &this_storage.Detector_3D.Array_p[loop_index][0][0], &this_storage.Detector_3D.Array_p2[loop_index][0][0], 778 part_filename); 779} 780 781%} 782 783FINALLY 784%{ 785// Remember to clean up allocated lists 786if (this_storage.temp_3DS_data.allocated_elements>0) free(this_storage.temp_3DS_data.elements); 787 788free3Ddouble_3DS(this_storage.Detector_3D.Array_N,n3,storage_N); 789free3Ddouble_3DS(this_storage.Detector_3D.Array_p,n3,storage_p); 790free3Ddouble_3DS(this_storage.Detector_3D.Array_p2,n3,storage_2p); 791 792if (accepted_processes.num_elements > 0) free(accepted_processes.elements); 793if (accepted_volumes.num_elements > 0) free(accepted_volumes.elements); 794 795%} 796 797END 798 799