     }
     printf("\n");
 
+    ////////////////////////////////////////////////////////////////////////
+    // Set the CUDA device. This is useful if your machine has multiple GPUs
+    // on it. 
+    ////////////////////////////////////////////////////////////////////////
+    cudaError_t error_id = cudaSetDevice(cuda_device);
+    if (error_id == cudaSuccess) {
+        printf( "cudaSetDevice returned %d\n-> %s\n", (int)error_id, cudaGetErrorString(error_id) );
+    }
+    else{
+        printf( "cudaSetDevice failed on Device %d!\n\n",cuda_device);
+        exit(-1);
+    }
+
     if(radec_input==1) int success = doCalcRaDec(infile0, infile1, outfile, silent_on_GPU_testing, scale_factor, nbins, hist_lower_range, hist_upper_range, hist_bin_width, log_binning_flag, two_different_files, conv_factor_angle);
 
     else  int success = doCalcMpc(infile0, infile1, outfile, silent_on_GPU_testing, scale_factor, nbins, hist_lower_range, hist_upper_range, hist_bin_width, log_binning_flag, two_different_files, conv_factor_angle);

    cudaFree(d_delta0);  
    cudaFree(d_alpha1);
    cudaFree(d_delta1);  
xs    cudaFree(dev_hist);

    return 0;
}  
//////////////////////////////////////////////////////////////////////


xs




    float conv_factor_angle = 57.2957795; // 180/pi // For if we need to convert arcdistance to arcsec or arcmin
    int radec_input = 1; // are we using ra/dec coords, or x/y/z coords? 
    bool silent_on_GPU_testing = false;
    int cuda_device = 0;

    ////////////////////////////////////////////////////////////////////////////////
    ////////////////////////////////////////////////////////////////////////////////

    while ((c = getopt(argc, argv, "ao:L:l:w:smSd:")) != -1) {
        switch(c) {
            case 'L':
                printf("L is set\n");
