/////////////////////////////////////////////////////////////////////////////////////////////////////// // // // This program "flip.c" is used to show improved reconstruction abilities of the triangulation // // code. It generates more precise reconstructed positions by compensating for a currently // // unresolved systematic error in the simulation code. This program compares the source position // // and reconstructed position of events to see if they are on the same side of the // // Ligo-Virgo-Hanford interferometer (IFO) plane. If not, the reconstructed position // // is flipped across the IFO plane. // // The program requires two inputs, the abbreviation of the waveform (SGQ3, SGQ9, CSGQ9, WNB) and // // the number of lines in the simulation results file to process. // // // /////////////////////////////////////////////////////////////////////////////////////////////////////// #include // Required c library files #include #include #include #define N1 0.2875 // The components of the IFO plane's normal vector. This is used to calculate which side of the IFO plane each event's source and recovered position is on. #define N2 -0.3704 #define N3 0.8833 // The vector_norm function takes three vector components and calculates that vector's norm. float vector_norm(float element_1, float element_2, float element_3) { float norm; norm = sqrt( (element_1*element_1) + (element_2*element_2) + (element_3*element_3)); return norm; } int main(int argc, char *argv[]) { if (argc == 3) // Check to see that there are two arguments to this function. { int number_of_events = atoi(argv[2]); FILE *data_file; // The file with output from our algorithm (recovered position, amplitude, energy). FILE *output_file; // The resulting file created for the prcfom.m script. char data_filename[40]; char output_filename[40]; /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // This section checks the name of files to be loaded and processed. This allows integration of this script into existing // shell scripts to automate the analysis of simulation results. if ( strcmp(argv[1], "SG") == 0 || strcmp(argv[1], "SGQ3") == 0 || strcmp(argv[1], "CSGQ9") == 0 || strcmp(argv[1], "WNB") == 0) // This checks the correct input file for retrieving data based on the input waveform. { sprintf(data_filename,"Final-PRC-%s.txt", argv[1]); data_file = fopen(data_filename,"r"); } else // If the filename does not contain a proper waveform abbreviation, notify the user and exit the program. { printf("\n\nYou must input a valid Waveform type (SG, SGQ3, CSGQ9, WNB)\n\n"); return 0; } /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Create a string with the name of the new file to be created, specifying that it holds the results of 'flip.c' sprintf(output_filename,"Final-PRC-%s-flipped.txt", argv[1]); // Open the 'flipped' results file for writing. output_file = fopen(output_filename,"w"); // Defining variables to store the data for each event (taken from the results file) char data_line[124]; char Injected_time[40]; char Amplitude[15]; char hrss[40]; int event_number; double Injected_Cos_Theta, Injected_Theta; float injected_theta, injected_phi, recovered_theta, recovered_phi, v1, v2, v3, r1, r2, r3, newr1, newr2, newr3; char status; // This for loop goes through the simulation results file line by line, flipping the reconstructed sources which are on opposite // sides of the IFO plane from their true source. The loop will end after looping through the specified number of lines in, or // reaching the end of, the simulation results file. for (event_number = 1; (event_number <= number_of_events) && (status != EOF) ; event_number++) { fgets(data_line,126,data_file); // This gets the next line of text from the data file, containing the results for 1 event. status = sscanf(data_line," %s %s %f %f %f %f %s \n", Injected_time, Amplitude, &injected_theta, &injected_phi, &recovered_theta, &recovered_phi, hrss); // This checks to see if the end of the data file has been reached. v1 = sin(injected_theta)*cos(injected_phi); // The components of the position vector of the event's source position v2 = sin(injected_theta)*sin(injected_phi); v3 = cos(injected_theta); r1 = sin(recovered_theta)*cos(recovered_phi); // The components of the position vector of the event's position recovered by the triangulation algorithm. r2 = sin(recovered_theta)*sin(recovered_phi); r3 = cos(recovered_theta); if ( ((v1*N1) + (v2*N2) + (v3*N3)) < 0 ) // If normal-dot-injected_position < 0 (on opposite sides of IFO plane) { if ( ((r1*N1) + (r2*N2) + (r3*N3)) > 0 ) // If normal-dot-recovered_position > 0 (on opp side of injected position, flip it) { newr1 = r1 - ( (2)*((r1*N1)+(r2*N2)+(r3*N3))*N1 ); // Calculate the components of the new recovered source position vector. newr2 = r2 - ( (2)*((r1*N1)+(r2*N2)+(r3*N3))*N2 ); newr3 = r3 - ( (2)*((r1*N1)+(r2*N2)+(r3*N3))*N3 ); recovered_phi = atan2(newr2,newr1); // These 2 lines calculate the flipped position and overwrite the original value. recovered_theta = acos(newr3/vector_norm(newr1, newr2, newr3)); } } else if ( ((v1*N1) + (v2*N2) + (v3*N3)) > 0 ) // Else, check to see if normal-dot-injected_position > 0 and normal-dot-recovered_position < 0. { if ( ((r1*N1) + (r2*N2) + (r3*N3)) < 0 ) { newr1 = r1 - ( (2)*((r1*N1)+(r2*N2)+(r3*N3))*N1 );; // Calculate the components of the new recovered source position vector. newr2 = r2 - ( (2)*((r1*N1)+(r2*N2)+(r3*N3))*N2 ); newr3 = r3 - ( (2)*((r1*N1)+(r2*N2)+(r3*N3))*N3 ); recovered_phi = atan2(newr2,newr1); // These 2 lines calculate the flipped position and overwrite the original value. recovered_theta = acos(newr3/vector_norm(newr1, newr2, newr3)); } } // If normal-dot-injected_position and normal-dot-recovered_position are of equal sign, then no change will take place. // Then, an new line will be written for each event, in the same format that is used in the original result file. fprintf(output_file,"%s\t%s\t%f\t%f\t%f\t%1.6f\t%s\t2\t2\n", Injected_time, Amplitude, injected_theta, injected_phi, recovered_theta, recovered_phi, hrss); } fclose(data_file); // Close the input and output files this program uses. fclose(output_file); } // If the number of program inputs does not equal two, print out a brief description of the inputs and exit the program. else printf("\nPlease input the waveform (SGQ3, SGQ9, CSGQ9, WNB), and the number of lines in the data file.\n"); return 0; }