/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright(C) 2007 Risoe National Laboratory. * * %I * Written by: Mads Bertelsen * Date: 20.08.15 * Version: $Revision: 0.1 $ * Origin: University of Copenhagen * * Three dimensional logger for spatial dimensions * * %D * Part of the Union components, a set of components that work together and thus * sperates geometry and physics within McStas. * The use of this component requires other components to be used. * * 1) One specifies a number of processes using process components * 2) These are gathered into material definitions using Union_make_material * 3) Geometries are placed using Union_box/cylinder/sphere, assigned a material * 4) A Union_master component placed after all of the above * * Only in step 4 will any simulation happen, and per default all geometries * defined before this master, but after the previous will be simulated here. * * There is a dedicated manual available for the Union_components * * This logger logs a 2D projection of the position of each scattering in the lab * frame. It does this in n3 space slices, so one can get a full 3D histogram of * scattering positions. * * A logger will log something for scattering events happening to certain volumes, * which are specified in the target_geometry string. By leaving it blank, all * geometries are logged, even the ones not defined at this point in the * instrument file. If a list og target_geometries is selected, one can further * narrow the events logged by providing a list of process names in target_process * which need to correspond with names of defined Union_process components. * * To use the logger_conditional_extend function, set it to some integer value n * and make and extend section to the master component that runs the geometry. * In this extend function, logger_conditional_extend[n] is 1 if the conditional * stack evaluated to true, 0 if not. This way one can check what rays is logged * using regular McStas monitors. Only works if a conditional is applied to this * logger. * * * %P * INPUT PARAMETERS: * D_direction_1: [string] Direction for first axis ("x", "y" or "z") * D1_max: [m] histogram boundery, max position value for first axis * D1_min: [m] histogram boundery, min position value for first axis * n1: [1] number of bins for first axis * D_direction_2: [string] Direction for second axis ("x", "y" or "z") * D2_max: [m] histogram boundery, max position value for second axis * D2_min: [m] histogram boundery, min position value for second axis * n2: [1] number of bins for second axis * D_direction_3: [string] Direction for second axis ("x", "y" or "z") * D3_max: [m] histogram boundery, max position value for third axis * D3_min: [m] histogram boundery, min position value for third axis * n3: [1] number of bins for third axis * filename: [string] Filename of produced data file * target_geometry: [string] Comma seperated list of geometry names that will be logged, leave empty for all volumes (even not defined yet) * target_process: [string] Comma seperated names of physical processes, if volumes are selected, one can select Union_process names * order_total: [1] Only log rays that scatter for the n'th time, 0 for all orders * order_volume: [1] Only log rays that scatter for the n'th time in the same geometry * order_volume_process: [1] Only log rays that scatter for the n'th time in the same geometry, using the same process * 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 * init: [string] name of Union_init component (typically "init", default) * * CALCULATED PARAMETERS: * * GLOBAL PARAMETERS: * * %L * * %E ******************************************************************************/ DEFINE COMPONENT Union_logger_3D_space SETTING 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", n1=90, n2=90, n3=10, order_total=0,order_volume=0,order_volume_process=0,logger_conditional_extend_index=-1, string init="init") /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ SHARE %{ #ifndef Union #error "The Union_init component must be included before this Union_logger_3D_space component" #endif struct temp_3DS_data_element_struct { int index_1; int index_2; int index_3; double weight; }; struct temp_3DS_data_struct { int num_elements; int allocated_elements; struct temp_3DS_data_element_struct* elements; }; struct a_3DS_storage_struct { struct Detector_3D_struct Detector_3D; struct temp_3DS_data_struct temp_3DS_data; int dim_1_choice; int dim_2_choice; int dim_3_choice; int order; int order_in_this_volume; int order_process_in_this_volume; Coords position; Rotation rotation; Rotation t_rotation; }; // record_to_temp // Would be nice if x y z, k_new and k_old were all coords void 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) { struct a_3DS_storage_struct* storage; storage = logger->data_union.p_3DS_storage; int add_point = 1; if (storage->order != 0) { if (storage->order - 1 == total_number_of_scattering_events) add_point = 1; else add_point = 0; } if (storage->order_in_this_volume != 0) { if (storage->order_in_this_volume - 1 == scattered_in_this_volume) add_point = 1; else add_point = 0; } if (storage->order_process_in_this_volume != 0) { if (storage->order_process_in_this_volume - 1 == scattered_in_this_volume_by_this_process) add_point = 1; else add_point = 0; } if (add_point == 1) { double p1, p2, p3; p1 = p2 = p3 = 0; // dim_1_choice = 0 for x, 1 for y, 2 for z. Done in initialize from input. "x" "y" "z". if (storage->dim_1_choice == 0) p1 = position->x - storage->position.x; else if (storage->dim_1_choice == 1) p1 = position->y - storage->position.y; else if (storage->dim_1_choice == 2) p1 = position->z - storage->position.z; if (storage->dim_2_choice == 0) p2 = position->x - storage->position.x; else if (storage->dim_2_choice == 1) p2 = position->y - storage->position.y; else if (storage->dim_2_choice == 2) p2 = position->z - storage->position.z; if (storage->dim_3_choice == 0) p2 = position->x - storage->position.x; else if (storage->dim_3_choice == 1) p2 = position->y - storage->position.y; else if (storage->dim_3_choice == 2) p2 = position->z - storage->position.z; int i, j, k; // Find bin in histogram 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) { i = floor ((p1 - storage->Detector_3D.D1min) * storage->Detector_3D.bins_1 / (storage->Detector_3D.D1max - storage->Detector_3D.D1min)); j = floor ((p2 - storage->Detector_3D.D2min) * storage->Detector_3D.bins_2 / (storage->Detector_3D.D2max - storage->Detector_3D.D2min)); k = floor ((p3 - storage->Detector_3D.D3min) * storage->Detector_3D.bins_3 / (storage->Detector_3D.D3max - storage->Detector_3D.D3min)); // Save bin in histogram to temp (may need to allocate more memory) int index; // 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); if (storage->temp_3DS_data.num_elements < storage->temp_3DS_data.allocated_elements) { storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_1 = k; storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_2 = i; storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_3 = j; storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements++].weight = p; } else { // No more space, need to allocate a larger buffer for this logger. Wish I had generics. // copy current data to temp struct temp_3DS_data_struct temporary_storage; temporary_storage.num_elements = storage->temp_3DS_data.num_elements; temporary_storage.elements = malloc (temporary_storage.num_elements * sizeof (struct temp_3DS_data_element_struct)); for (index = 0; index < storage->temp_3DS_data.num_elements; index++) { temporary_storage.elements[index].index_1 = storage->temp_3DS_data.elements[index].index_1; temporary_storage.elements[index].index_2 = storage->temp_3DS_data.elements[index].index_2; temporary_storage.elements[index].index_3 = storage->temp_3DS_data.elements[index].index_3; temporary_storage.elements[index].weight = storage->temp_3DS_data.elements[index].weight; } // free current data free (storage->temp_3DS_data.elements); // allocate larger array (10 larger) storage->temp_3DS_data.allocated_elements = 10 + storage->temp_3DS_data.num_elements; storage->temp_3DS_data.elements = malloc (storage->temp_3DS_data.allocated_elements * sizeof (struct temp_3DS_data_element_struct)); // copy back from temp for (index = 0; index < storage->temp_3DS_data.num_elements; index++) { storage->temp_3DS_data.elements[index].index_1 = temporary_storage.elements[index].index_1; storage->temp_3DS_data.elements[index].index_2 = temporary_storage.elements[index].index_2; storage->temp_3DS_data.elements[index].index_3 = temporary_storage.elements[index].index_3; storage->temp_3DS_data.elements[index].weight = temporary_storage.elements[index].weight; } // free temporary data free (temporary_storage.elements); // add new data point storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_1 = k; storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_2 = i; storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements].index_3 = j; storage->temp_3DS_data.elements[storage->temp_3DS_data.num_elements++].weight = p; } // 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 if (storage->temp_3DS_data.num_elements == 1) add_to_logger_with_data (logger_with_data_array, logger); } } } // clear_temp void clear_temp_3DS (union logger_data_union* data_union) { data_union->p_3DS_storage->temp_3DS_data.num_elements = 0; } // record_to_perm void 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) { // printf("In record to permanent \n"); struct a_3DS_storage_struct* storage; storage = logger->data_union.p_3DS_storage; int add_point = 1; if (storage->order != 0) { if (storage->order - 1 == total_number_of_scattering_events) add_point = 1; else add_point = 0; } if (storage->order_in_this_volume != 0) { if (storage->order_in_this_volume - 1 == scattered_in_this_volume) add_point = 1; else add_point = 0; } if (storage->order_process_in_this_volume != 0) { if (storage->order_process_in_this_volume - 1 == scattered_in_this_volume_by_this_process) add_point = 1; else add_point = 0; } if (add_point == 1) { // printf("storage was set \n"); double p1, p2, p3; p1 = p2 = p3 = 0; // dim_1_choice = 0 for x, 1 for y, 2 for z. Done in initialize from input. "x" "y" "z". if (storage->dim_1_choice == 0) p1 = position->x - storage->position.x; else if (storage->dim_1_choice == 1) p1 = position->y - storage->position.y; else if (storage->dim_1_choice == 2) p1 = position->z - storage->position.z; if (storage->dim_2_choice == 0) p2 = position->x - storage->position.x; else if (storage->dim_2_choice == 1) p2 = position->y - storage->position.y; else if (storage->dim_2_choice == 2) p2 = position->z - storage->position.z; if (storage->dim_3_choice == 0) p3 = position->x - storage->position.x; else if (storage->dim_3_choice == 1) p3 = position->y - storage->position.y; else if (storage->dim_3_choice == 2) p3 = position->z - storage->position.z; int i, j, k; // Find bin in histogram 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) { i = floor ((p1 - storage->Detector_3D.D1min) * storage->Detector_3D.bins_1 / (storage->Detector_3D.D1max - storage->Detector_3D.D1min)); j = floor ((p2 - storage->Detector_3D.D2min) * storage->Detector_3D.bins_2 / (storage->Detector_3D.D2max - storage->Detector_3D.D2min)); k = floor ((p3 - storage->Detector_3D.D3min) * storage->Detector_3D.bins_3 / (storage->Detector_3D.D3max - storage->Detector_3D.D3min)); /* printf("(p1,p2,p3) = (%f,%f,%f)\n",p1,p2,p3); printf("limits1 = [%f,%f] \n",storage->Detector_3D.D1min,storage->Detector_3D.D1max); printf("limits2 = [%f,%f] \n",storage->Detector_3D.D2min,storage->Detector_3D.D2max); printf("limits3 = [%f,%f] \n",storage->Detector_3D.D3min,storage->Detector_3D.D3max); 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)); printf("Added to statistics for monitor [%d] [%d] [%d] \n",i,j,k); printf("indicies found\n"); */ /* storage->Detector_3D.Array_N[i][j][k]++; storage->Detector_3D.Array_p[i][j][k] += p; storage->Detector_3D.Array_p2[i][j][k] += p*p; */ // because of the order in which the elements are laid out in memory, the k index must be first. storage->Detector_3D.Array_N[k][i][j]++; storage->Detector_3D.Array_p[k][i][j] += p; storage->Detector_3D.Array_p2[k][i][j] += p * p; } } } // write_temp_to_perm void write_temp_to_perm_3DS (union logger_data_union* data_union) { struct a_3DS_storage_struct* storage; storage = data_union->p_3DS_storage; int index; // Add all data points to the historgram, they are saved as index / weight combinations for (index = 0; index < storage->temp_3DS_data.num_elements; index++) { 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]++; 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; 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; } clear_temp_3DS (data_union); } void write_temp_to_perm_final_p_3DS (union logger_data_union* data_union, double final_weight) { struct a_3DS_storage_struct* storage; storage = data_union->p_3DS_storage; int index; // Add all data points to the historgram, they are saved as index / weight combinations for (index = 0; index < storage->temp_3DS_data.num_elements; index++) { 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]++; 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; 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; } clear_temp_3DS (data_union); } // Only need to define linking function for loggers once. #ifndef UNION_LOGGER #define UNION_LOGGER Dummy // Linking function for loggers, finds the indicies of the specified geometries on the global_geometry_list void 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) { // Need to check a input_string of text for an occurance of name. If it is in the inputstring, yes return 1, otherwise 0. char* token; int loop_index; char local_string[512]; strcpy (local_string, input_string); // get the first token token = strtok (local_string, ","); // walk through other tokens while (token != NULL) { // printf( " %s\n", token ); for (loop_index = 0; loop_index < global_geometry_list->num_elements; loop_index++) { if (strcmp (token, global_geometry_list->elements[loop_index].name) == 0) { add_element_to_int_list (accepted_volumes, loop_index); break; } if (loop_index == global_geometry_list->num_elements - 1) { // All possible geometry names have been looked through, and the break was not executed. // Alert the user to this problem by showing the geometry name that was not found and the currently available geometires printf ("\n"); 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); printf (" The unrecoignized geometry name was: \"%s\" \n", token); printf (" The geometries available at this point (need to be defined before the logger): \n"); for (loop_index = 0; loop_index < global_geometry_list->num_elements; loop_index++) printf (" %s\n", global_geometry_list->elements[loop_index].name); exit (1); } } // Updates the token token = strtok (NULL, ","); } } void 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) { // Need to check a input_string of text for an occurance of name. If it is in the inputstring, yes return 1, otherwise 0. char* token; int loop_index; char local_string[256]; strcpy (local_string, input_string); // get the first token token = strtok (local_string, ","); // walk through other tokens while (token != NULL) { // printf( " %s\n", token ); for (loop_index = 0; loop_index < p_physics->number_of_processes; loop_index++) { if (strcmp (token, p_physics->p_scattering_array[loop_index].name) == 0) { add_element_to_int_list (accepted_processes, loop_index); break; } if (loop_index == p_physics->number_of_processes - 1) { // All possible process names have been looked through, and the break was not executed. // Alert the user to this problem by showing the process name that was not found and the currently available processes printf ("\n"); 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); printf (" The unrecoignized process name was: \"%s\" \n", token); printf (" The processes defined in material \"%s\" of which Volume \"%s\" is made: \n", p_physics->name, Volume_name); for (loop_index = 0; loop_index < p_physics->number_of_processes; loop_index++) printf (" %s\n", p_physics->p_scattering_array[loop_index].name); exit (1); } } // Updates the token token = strtok (NULL, ","); } } #endif double*** allocate3Ddouble_3DS (int count_x, int count_y, int count_z, double* storage) { // double *storage = malloc(count_x * count_y * count_z * sizeof(double)); storage = malloc (count_x * count_y * count_z * sizeof (double)); if (storage == NULL) { free (storage); printf ("\nERROR: storage array not allocated in Union_logger_3D_space\n"); } double* alloc = storage; double*** x; int i, j; x = malloc (count_x * sizeof (*x)); if (x == NULL) { free (x); printf ("\nERROR: x array not allocated in Union_logger_3D_space\n"); exit (1); } for (i = 0; i < count_x; i++) { x[i] = malloc (count_y * sizeof (**x)); if (x[i] == NULL) { free (x[i]); printf ("\nERROR: x[i] array not allocated in Union_logger_3D_space\n"); exit (1); } for (j = 0; j < count_y; j++) { x[i][j] = alloc; alloc += count_z; } } return x; } void free3Ddouble_3DS (double*** ptr_array, int count_x, double* storage) { if (!ptr_array) return; int x; free (storage); for (x = 0; x < count_x; x++) { free (ptr_array[x]); } free (ptr_array); } %} DECLARE %{ // From make material // Needed for transport to the main component // struct global_material_element_struct global_material_element; // struct physics_struct this_material; int loop_index; int process_index; int specified_processes; char local_string[256]; char number_string[16]; char part_filename[256]; // Reused for logger struct pointer_to_1d_int_list accepted_processes; struct global_logger_element_struct logger_list_element; struct pointer_to_1d_int_list accepted_volumes; struct logger_struct this_logger; struct a_3DS_storage_struct this_storage; struct loggers_struct* loggers_on_target_volume; struct Volume_struct* target_volume; double* storage_N; double* storage_p; double* storage_2p; %} INITIALIZE %{ accepted_processes.num_elements = 0; accepted_processes.elements = NULL; accepted_volumes.num_elements = 0; accepted_volumes.elements = NULL; // Initialize storage from input if (D1_min >= D1_max) { printf ("ERROR, Union logger \"%s\" had D1_min >= D1_max.\n", NAME_CURRENT_COMP); exit (1); } this_storage.Detector_3D.D1min = D1_min; this_storage.Detector_3D.D1max = D1_max; if (D2_min >= D2_max) { printf ("ERROR, Union logger \"%s\" had D2_min >= D2_max.\n", NAME_CURRENT_COMP); exit (1); } this_storage.Detector_3D.D2min = D2_min; this_storage.Detector_3D.D2max = D2_max; if (D3_min >= D3_max) { printf ("ERROR, Union logger \"%s\" had D3_min >= D3_max.\n", NAME_CURRENT_COMP); exit (1); } this_storage.Detector_3D.D3min = D3_min; this_storage.Detector_3D.D3max = D3_max; if (n1 <= 0) { printf ("ERROR, Union logger \"%s\" had n1 <= 0.\n", NAME_CURRENT_COMP); exit (1); } this_storage.Detector_3D.bins_1 = n1; if (n2 <= 0) { printf ("ERROR, Union logger \"%s\" had n2 <= 0.\n", NAME_CURRENT_COMP); exit (1); } this_storage.Detector_3D.bins_2 = n2; if (n3 <= 0) { printf ("ERROR, Union logger \"%s\" had n3 <= 0.\n", NAME_CURRENT_COMP); exit (1); } this_storage.Detector_3D.bins_3 = n3; // printf("past input sanitation \n"); // Remember to take special care when deallocating this array, use free3Ddouble /* this_storage.Detector_3D.Array_N = allocate3Ddouble_3DS(n1,n2,n3); // Here the n1 double is cast to an int this_storage.Detector_3D.Array_p = allocate3Ddouble_3DS(n1,n2,n3); this_storage.Detector_3D.Array_p2 = allocate3Ddouble_3DS(n1,n2,n3); */ // D3 is in poisition 1, because that gives continous memory in XY plane for easy plotting this_storage.Detector_3D.Array_N = allocate3Ddouble_3DS (n3, n1, n2, storage_N); // Here the n1 double is cast to an int this_storage.Detector_3D.Array_p = allocate3Ddouble_3DS (n3, n1, n2, storage_p); this_storage.Detector_3D.Array_p2 = allocate3Ddouble_3DS (n3, n1, n2, storage_2p); /* // Error in the order? this_storage.Detector_3D.Array_N = allocate3Ddouble_3DS(n1,n3,n2); // Here the n1 double is cast to an int this_storage.Detector_3D.Array_p = allocate3Ddouble_3DS(n1,n3,n2); this_storage.Detector_3D.Array_p2 = allocate3Ddouble_3DS(n1,n3,n2); */ // printf("Allocated 3D arrays \n"); int l1, l2, l3; for (l1 = 0; l1 < n1; l1++) { // n1 is technically a double, but this works fine for (l2 = 0; l2 < n2; l2++) { for (l3 = 0; l3 < n3; l3++) { this_storage.Detector_3D.Array_N[l3][l1][l2] = 0; this_storage.Detector_3D.Array_p[l3][l1][l2] = 0; this_storage.Detector_3D.Array_p2[l3][l1][l2] = 0; } } } // printf("Initialized 3D arrays \n"); // printf("past 3D pointer assignment \n"); // Input sanitation for filename apparently done in 3D_detector_out if (strcmp (D_direction_1, "x") == 0 || strcmp (D_direction_1, "X") == 0) { this_storage.dim_1_choice = 0; sprintf (this_storage.Detector_3D.string_axis_1, "x [m]"); sprintf (this_storage.Detector_3D.title_string, "3D position Union logger in plane x"); } else if (strcmp (D_direction_1, "y") == 0 || strcmp (D_direction_1, "Y") == 0) { this_storage.dim_1_choice = 1; sprintf (this_storage.Detector_3D.string_axis_1, "y [m]"); sprintf (this_storage.Detector_3D.title_string, "3D position Union logger in plane y"); } else if (strcmp (D_direction_1, "z") == 0 || strcmp (D_direction_1, "Z") == 0) { this_storage.dim_1_choice = 2; sprintf (this_storage.Detector_3D.string_axis_1, "z [m]"); sprintf (this_storage.Detector_3D.title_string, "3D position Union logger in plane z"); } else { 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); exit (1); } char temp_string[2] = "UN"; // UN stands for unset, courtesy of linter play if (strcmp (D_direction_2, "x") == 0 || strcmp (D_direction_2, "X") == 0) { this_storage.dim_2_choice = 0; sprintf (this_storage.Detector_3D.string_axis_2, "x [m]"); sprintf (temp_string, "x"); } else if (strcmp (D_direction_2, "y") == 0 || strcmp (D_direction_2, "Y") == 0) { this_storage.dim_2_choice = 1; sprintf (this_storage.Detector_3D.string_axis_2, "y [m]"); sprintf (temp_string, "y"); } else if (strcmp (D_direction_2, "z") == 0 || strcmp (D_direction_2, "Z") == 0) { this_storage.dim_2_choice = 2; sprintf (this_storage.Detector_3D.string_axis_2, "z [m]"); sprintf (temp_string, "z"); } else { 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); exit (1); } 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 if (strcmp (D_direction_3, "x") == 0 || strcmp (D_direction_3, "X") == 0) { this_storage.dim_3_choice = 0; sprintf (this_storage.Detector_3D.string_axis_3, "x [m]"); sprintf (temp_string, "x"); } else if (strcmp (D_direction_3, "y") == 0 || strcmp (D_direction_3, "Y") == 0) { this_storage.dim_3_choice = 1; sprintf (this_storage.Detector_3D.string_axis_3, "y [m]"); sprintf (temp_string, "y"); } else if (strcmp (D_direction_3, "z") == 0 || strcmp (D_direction_3, "Z") == 0) { this_storage.dim_3_choice = 2; sprintf (this_storage.Detector_3D.string_axis_3, "z [m]"); sprintf (temp_string, "z"); } else { 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); exit (1); } strcat (this_storage.Detector_3D.title_string, temp_string); sprintf (this_storage.Detector_3D.Filename, "%s", filename); this_storage.order = order_total; this_storage.order_in_this_volume = order_volume; this_storage.order_process_in_this_volume = order_volume_process; this_storage.temp_3DS_data.num_elements = 0; this_storage.temp_3DS_data.allocated_elements = 10; this_storage.temp_3DS_data.elements = malloc (this_storage.temp_3DS_data.allocated_elements * sizeof (struct temp_3DS_data_element_struct)); if (_getcomp_index (init) < 0) { fprintf (stderr, "Union_logger_3D:%s: Error identifying Union_init component, %s is not a known component name.\n", NAME_CURRENT_COMP, init); exit (-1); } // Global variables retrieved from init component struct global_positions_to_transform_list_struct* global_positions_to_transform_list = COMP_GETPAR3 (Union_init, init, global_positions_to_transform_list); struct global_rotations_to_transform_list_struct* global_rotations_to_transform_list = COMP_GETPAR3 (Union_init, init, global_rotations_to_transform_list); this_storage.position = POS_A_CURRENT_COMP; add_position_pointer_to_list (global_positions_to_transform_list, &this_storage.position); rot_copy (this_storage.rotation, ROT_A_CURRENT_COMP); add_rotation_pointer_to_list (global_rotations_to_transform_list, &this_storage.rotation); rot_transpose (ROT_A_CURRENT_COMP, this_storage.t_rotation); add_rotation_pointer_to_list (global_rotations_to_transform_list, &this_storage.t_rotation); // printf("past direction choices sanitation \n"); // Book keeping this_logger.logger_extend_index = logger_conditional_extend_index; this_logger.function_pointers.active_record_function = &record_to_perm_3DS; // Assume no conditional this_logger.function_pointers.inactive_record_function = &record_to_temp_3DS; // If an assume is present, these two pointers are switched // Temp to perm functions, and standard identifier this_logger.function_pointers.select_t_to_p = 1; // 1: temp_to_perm, 2: temp_to_perm_final_p this_logger.function_pointers.temp_to_perm = &write_temp_to_perm_3DS; this_logger.function_pointers.temp_to_perm_final_p = &write_temp_to_perm_final_p_3DS; this_logger.function_pointers.clear_temp = &clear_temp_3DS; // Initializing for conditional this_logger.conditional_list.num_elements = 0; // this_logger.function_pointers.perm_to_disk = &write_perm_to_disk_3DS; //Obsolete // printf("past this_logger function assignment \n"); this_logger.data_union.p_3DS_storage = &this_storage; sprintf (this_logger.name, "%s", NAME_CURRENT_COMP); // printf("past this_logger assignment \n"); sprintf (logger_list_element.name, "%s", NAME_CURRENT_COMP); logger_list_element.component_index = INDEX_CURRENT_COMP; logger_list_element.logger = &this_logger; // printf("past logger_list_element assignment \n"); // Get global variables from init component struct pointer_to_global_geometry_list* global_geometry_list = COMP_GETPAR3 (Union_init, init, global_geometry_list); struct pointer_to_global_logger_list* global_specific_volumes_logger_list = COMP_GETPAR3 (Union_init, init, global_specific_volumes_logger_list); // In order to run the logger at the right times, pointers to this logger is stored in each volume it logs, // and additionally for each avaiable process. If a process is not logged, the pointer is simply not stored. // Need to find the volumes for which the processes should have a reference to this logger if (target_geometry && strlen (target_geometry) && strcmp (target_geometry, "NULL") && strcmp (target_geometry, "0")) { // Certain volumes were selected, find the indicies in the global_geometry_list manual_linking_function_logger_volumes (target_geometry, global_geometry_list, &accepted_volumes, NAME_CURRENT_COMP); // Add this logger to the global_specific_volumes_logger_list (so that conditionals can affect it) add_element_to_logger_list (global_specific_volumes_logger_list, logger_list_element); for (loop_index = 0; loop_index < accepted_volumes.num_elements; loop_index++) { target_volume = global_geometry_list->elements[accepted_volumes.elements[loop_index]].Volume; // Add an element to its logger list add_initialized_logger_in_volume (&target_volume->loggers, target_volume->p_physics->number_of_processes); if (target_process && strlen (target_process) && strcmp (target_process, "NULL") && strcmp (target_process, "0")) { // Unused process pointers should point to NULL for (process_index = 0; process_index < target_volume->p_physics->number_of_processes; process_index++) { target_volume->loggers.p_logger_volume[target_volume->loggers.num_elements - 1].p_logger_process[process_index] = NULL; } // A target_process was set, find it within the volume structure (can be many processes) manual_linking_function_logger_processes (target_process, target_volume->p_physics, &accepted_processes, NAME_CURRENT_COMP, target_volume->name); for (process_index = 0; process_index < accepted_processes.num_elements; process_index++) { // Add pointer to this logger for all the accepted processes in this newly added loggers element target_volume->loggers.p_logger_volume[target_volume->loggers.num_elements - 1].p_logger_process[accepted_processes.elements[process_index]] = &this_logger; } } else { // No target_process was set, attatch the logger to all processes for (process_index = 0; process_index < target_volume->p_physics->number_of_processes; process_index++) { target_volume->loggers.p_logger_volume[target_volume->loggers.num_elements - 1].p_logger_process[process_index] = &this_logger; } } } } else { // Send to global_all_volumes_logger_list // Here there is no system for selecting processes as well struct pointer_to_global_logger_list* global_all_volume_logger_list = COMP_GETPAR3 (Union_init, init, global_all_volume_logger_list); add_element_to_logger_list (global_all_volume_logger_list, logger_list_element); } %} TRACE %{ %} SAVE %{ // Write to disk for (loop_index = 0; loop_index < this_storage.Detector_3D.bins_3; loop_index++) { sprintf (number_string, "%d", loop_index); sprintf (part_filename, "%s_%s", this_storage.Detector_3D.Filename, number_string); DETECTOR_OUT_2D (this_storage.Detector_3D.title_string, this_storage.Detector_3D.string_axis_1, this_storage.Detector_3D.string_axis_2, this_storage.Detector_3D.D1min, this_storage.Detector_3D.D1max, this_storage.Detector_3D.D2min, this_storage.Detector_3D.D2max, this_storage.Detector_3D.bins_1, this_storage.Detector_3D.bins_2, &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], part_filename); } %} FINALLY %{ // Remember to clean up allocated lists if (this_storage.temp_3DS_data.allocated_elements > 0) free (this_storage.temp_3DS_data.elements); free3Ddouble_3DS (this_storage.Detector_3D.Array_N, n3, storage_N); free3Ddouble_3DS (this_storage.Detector_3D.Array_p, n3, storage_p); free3Ddouble_3DS (this_storage.Detector_3D.Array_p2, n3, storage_2p); if (accepted_processes.num_elements > 0) free (accepted_processes.elements); if (accepted_volumes.num_elements > 0) free (accepted_volumes.elements); %} END