- Timestamp:
- 09/23/03 17:50:55 (21 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
r2351 r2352 21 21 // 22 22 // $RCSfile: camera.cxx,v $ 23 // $Revision: 1. 59$23 // $Revision: 1.60 $ 24 24 // $Author: blanch $ 25 // $Date: 2003-09-2 2 11:38:47$25 // $Date: 2003-09-23 16:50:55 $ 26 26 // 27 27 //////////////////////////////////////////////////////////////////////// … … 134 134 @"*/ 135 135 136 /*!@" 137 138 And this is the information about the whole telescope. 139 140 @"*/ 141 136 142 //!@{ 137 static int ct_Type; //@< Type of telescope: 0:CT1, 1:MAGIC138 //!@}139 140 /*!@"141 142 And this is the information about the whole telescope.143 144 @"*/145 146 //!@{147 143 148 144 // parameters of the CT (from the CT definition file) 149 150 ////@: Focal distances [cm]151 //static float *ct_Focal;152 153 //@: Mean Focal distances [cm]154 static float ct_Focal_mean;155 156 //@: STDev. Focal distances [cm]157 static float ct_Focal_std;158 159 //@: Mean Point Spread function [cm]160 static float ct_PSpread_mean;161 162 //@: STDev. Point Spread function [cm]163 static float ct_PSpread_std;164 165 //@: STDev. Adjustmente deviation [cm]166 static float ct_Adjustment_std;167 168 //@: Radius of the Black Spot in mirror [cm]169 static float ct_BlackSpot_rad;170 171 //@: Radius of one mirror [cm]172 static float ct_RMirror;173 174 //@: Camera width [cm]175 static float ct_CameraWidth;176 177 //@: Pixel width [cm]178 static float ct_PixelWidth;179 180 //@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60)181 static float ct_PixelWidth_corner_2_corner;182 183 //@: Number of mirrors184 static int ct_NMirrors = 0;185 145 186 146 //@: Number of pixels 187 147 static int ct_NPixels; 188 189 //@: Number of pixels190 static int ct_NCentralPixels;191 192 //@: Number of pixels193 static int ct_NGapPixels;194 195 //@: name of the CT definition file to use196 static char ct_filename[256];197 148 198 149 //@: Number of CT … … 286 237 287 238 288 //@: coordinates x,y for each pixel289 static float **pixary;290 291 239 //@: contents of the pixels (ph.e.) 292 240 static float *fnpix; … … 331 279 332 280 @"*/ 333 334 //!@{335 // Pointer to a table with the following info.:336 337 static float **ct_data;338 339 /*340 * TYPE=0 (CT1)341 * i s rho theta x y z thetan phin xn yn zn342 *343 * i : number of the mirror344 * s : arc length [cm]345 * rho : polar rho of the position of the center of the mirror [cm]346 * theta : polar angle of the position of the center of the mirror [cm]347 * x : x coordinate of the center of the mirror [cm]348 * y : y coordinate of the center of the mirror [cm]349 * z : z coordinate of the center of the mirror [cm]350 * thetan : polar theta angle of the direction where the mirror points to351 * phin : polar phi angle of the direction where the mirror points to352 * xn : xn coordinate of the normal vector in the center (normalized)353 * yn : yn coordinate of the normal vector in the center (normalized)354 * zn : zn coordinate of the normal vector in the center (normalized)355 *356 * TYPE=1 (MAGIC)357 * i f sx sy x y z thetan phin358 *359 * i : number of the mirror360 * f : focal distance of that mirror361 * sx : curvilinear coordinate of mirror's center in X[cm]362 * sy : curvilinear coordinate of mirror's center in X[cm]363 * x : x coordinate of the center of the mirror [cm]364 * y : y coordinate of the center of the mirror [cm]365 * z : z coordinate of the center of the mirror [cm]366 * thetan : polar theta angle of the direction where the mirror points to367 * phin : polar phi angle of the direction where the mirror points to368 * xn : xn coordinate of the normal vector in the center (normalized)369 * yn : yn coordinate of the normal vector in the center (normalized)370 * zn : zn coordinate of the normal vector in the center (normalized)371 */372 //!@}373 281 374 282 /*!@" … … 737 645 strcpy( rootname, get_root_filename() ); 738 646 strcpy( rootname_loop, get_loop_filename() ); 739 strcpy( ct_filename, get_ct_filename() );740 647 strcpy( nsbpathname, get_nsb_directory() ); 741 648 strcpy( nsbpath_outer, get_nsb_directory_outer() ); … … 753 660 754 661 Starfield_rotate = get_starfield_rotate(); 755 log(SIGNATURE,756 "%s:\n\t%20s: %s\n",757 "Starfield Rotate",758 "Rotate Starfield", ONoff(Starfield_rotate)759 );760 662 761 663 // log filenames information 762 664 for(Int_t ict=0;ict<ct_Number;ict++){ 665 666 log(SIGNATURE,"\t%s : %i\n\t%20s:\t%i\n", 667 "------- TELESCOPE ",ict+1, 668 "Geometry ", GeometryCamera[ict]); 669 strcpy( qe_filename, get_qe_filename(ict)); 670 671 log(SIGNATURE, 672 "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n", 673 "Filenames", 674 "In", inname_CT[ict], 675 "Stars", starfieldname, 676 "NSB database","(inner pixels)", nsbpathname, 677 "(outer pixels)", nsbpath_outer, 678 "QE", qe_filename, 679 "Data", datname, 680 "ROOT", rootname 681 ); 682 683 // log Trigger information 684 685 if (Trigger_Loop) { 686 log(SIGNATURE, 687 "%s:\n\t%20s: from %5.2f to %5.2f by %5.2f step\n\t%20s: %i - %i\n\t%20s: %i - %i\n\t%20s\n", 688 "Trigger Loop mode", 689 "Threshold",Trigger_loop_lthres,Trigger_loop_uthres,Trigger_loop_sthres, 690 "Multiplicity",Trigger_loop_lmult,Trigger_loop_umult, 691 "Topology",Trigger_loop_ltop,Trigger_loop_utop, 692 rootname_loop); 693 } 694 else if (Individual_Thres_Pixel == FALSE){ 695 log(SIGNATURE, 696 "%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n", 697 "Single Trigger mode", 698 "Threshold",qThreshold[ict][0], 699 "Multiplicity",Trigger_multiplicity[ict], 700 "Topology",Trigger_topology[ict]); 701 } 702 else{ 703 log(SIGNATURE, 704 "%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n", 705 "Single Trigger mode", 706 "Threshold","Different for each pixel", 707 "Multiplicity",Trigger_multiplicity[ict], 708 "Topology",Trigger_topology[ict]); 709 } 710 log(SIGNATURE,"\t%s\n", 711 "END TELESCOPE --------"); 712 } 713 // log flags information 714 763 715 log(SIGNATURE, 764 "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n", 765 "Filenames", 766 "In", inname_CT[0], 767 "Stars", starfieldname, 768 "NSB database (inner pixels)", nsbpathname, 769 "NSB database (outer pixels)", nsbpath_outer, 770 "CT", ct_filename, 771 "Data", datname, 772 "ROOT", rootname 773 ); 774 775 // log Trigger information 776 777 if (Trigger_Loop) { 778 log(SIGNATURE, 779 "%s:\n\t%20s: from %5.2f to %5.2f by %5.2f step\n\t%20s: %i - %i\n\t%20s: %i - %i\n\t%20s\n", 780 "Trigger Loop mode", 781 "Threshold",Trigger_loop_lthres,Trigger_loop_uthres,Trigger_loop_sthres, 782 "Multiplicity",Trigger_loop_lmult,Trigger_loop_umult, 783 "Topology",Trigger_loop_ltop,Trigger_loop_utop, 784 rootname_loop); 785 } 786 else if (Individual_Thres_Pixel == FALSE){ 787 log(SIGNATURE, 788 "%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n", 789 "Single Trigger mode", 790 "Threshold",qThreshold[0][0], 791 "Multiplicity",Trigger_multiplicity[0], 792 "Topology",Trigger_topology[0]); 793 } 794 else{ 795 log(SIGNATURE, 796 "%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n", 797 "Single Trigger mode", 798 "Threshold","Different for each pixel", 799 "Multiplicity",Trigger_multiplicity[0], 800 "Topology",Trigger_topology[0]); 801 } 802 // log flags information 803 804 log(SIGNATURE, 805 "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n", 716 "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %3.2f(%s) %3.2f(%s) %s\n", 806 717 "Flags", 807 718 "Data_From_STDIN", ONoff(Data_From_STDIN), 808 719 "Write_All_Events", ONoff(Write_All_Images), 809 "Write_All_Data", ONoff(Write_All_Data)); 720 "Write_All_Data", ONoff(Write_All_Data), 721 "Rotate Starfield", ONoff(Starfield_rotate), 722 "Electronic Noise", Trigger_noise,"trigger", FADC_noise,"fadc",ONoff(addElecNoise) 723 ); 810 724 811 725 // log flags information … … 822 736 823 737 log(SIGNATURE, 824 "%s:\n\t%20s: %f \n",738 "%s:\n\t%20s: %f: %s\n\t%20s: %f\n", 825 739 "Parameters", 826 "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB)); 740 "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB), 741 "Pedestals = ", get_FADC_pedestal() 742 ); 827 743 828 744 // log selections … … 881 797 // read parameters from the ct.def file 882 798 883 read_ct_file();799 //read_ct_file(); 884 800 885 801 Int_t Lev0, Lev1; … … 2791 2707 2792 2708 //!----------------------------------------------------------- 2793 // @name read_ct_file2794 //2795 // @desc read CT definition file2796 //2797 // @date Sat Jun 27 05:58:56 MET DST 19982798 //------------------------------------------------------------2799 // @function2800 2801 //!@{2802 void2803 read_ct_file(void)2804 {2805 char line[LINE_MAX_LENGTH]; //@< line to get from the ctin2806 char token[ITEM_MAX_LENGTH]; //@< a single token2807 int i, j; //@< dummy counters2808 2809 log( "read_ct_file", "start.\n" );2810 2811 ifstream ctin ( ct_filename );2812 2813 if ( ctin.bad() )2814 error( "read_ct_file",2815 "Cannot open CT def. file: %s\n", ct_filename );2816 2817 // loop till the "end" directive is reached2818 2819 while (!ctin.eof()) {2820 2821 // get line from stdin2822 2823 ctin.getline(line, LINE_MAX_LENGTH);2824 2825 // look for each item at the beginning of the line2826 2827 for (i=0; i<=define_mirrors; i++)2828 if (strstr(line, CT_ITEM_NAMES[i]) == line)2829 break;2830 2831 // if it is not a valid line, just ignore it2832 2833 if (i == define_mirrors+1)2834 continue;2835 2836 // case block for each directive2837 2838 switch ( i ) {2839 2840 case type: // <type of telescope> (0:CT1 ¦ 1:MAGIC)2841 2842 // get focal distance2843 2844 sscanf(line, "%s %d", token, &ct_Type);2845 2846 log( "read_ct_file", "<Type of Telescope>: %s\n",2847 ((ct_Type==0) ? "CT1" : "MAGIC") );2848 2849 break;2850 2851 case focal_distance: // <focal distance> [cm]2852 2853 // get focal distance2854 2855 sscanf(line, "%s %f", token, &ct_Focal_mean);2856 2857 log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );2858 2859 break;2860 2861 case focal_std: // s(focal distance) [cm]2862 2863 // get focal distance2864 2865 sscanf(line, "%s %f", token, &ct_Focal_std);2866 2867 log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );2868 2869 break;2870 2871 case point_spread: // <point spread> [cm]2872 2873 // get point spread2874 2875 sscanf(line, "%s %f", token, &ct_PSpread_mean);2876 2877 log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );2878 2879 break;2880 2881 case point_std: // s(point spread) [cm]2882 2883 // get point spread2884 2885 sscanf(line, "%s %f", token, &ct_PSpread_std);2886 2887 log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );2888 2889 break;2890 2891 case adjustment_dev: // s(adjustment_dev) [cm]2892 2893 // get point spread2894 2895 sscanf(line, "%s %f", token, &ct_Adjustment_std);2896 2897 log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );2898 2899 break;2900 2901 case black_spot: // radius of the black spot in the center [cm]2902 2903 // get black spot radius2904 2905 sscanf(line, "%s %f", token, &ct_BlackSpot_rad);2906 2907 log( "read_ct_file", "Radius of the black spots: %f cm\n",2908 ct_BlackSpot_rad);2909 2910 break;2911 2912 case r_mirror: // radius of the mirrors [cm]2913 2914 // get radius of mirror2915 2916 sscanf(line, "%s %f", token, &ct_RMirror);2917 2918 log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );2919 2920 break;2921 2922 case n_mirrors: // number of mirrors2923 2924 // get the name of the output_file from the line2925 2926 sscanf(line, "%s %d", token, &ct_NMirrors);2927 2928 log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );2929 2930 break;2931 2932 case camera_width: // camera width [cm]2933 2934 // get the name of the ct_file from the line2935 2936 sscanf(line, "%s %f", token, &ct_CameraWidth);2937 2938 log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );2939 2940 break;2941 2942 case n_pixels: // number of pixels2943 2944 // get the name of the output_file from the line2945 2946 sscanf(line, "%s %d", token, &ct_NPixels);2947 2948 log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );2949 2950 break;2951 2952 case n_centralpixels: // number of central pixels2953 2954 // get the name of the output_file from the line2955 2956 sscanf(line, "%s %d", token, &ct_NCentralPixels);2957 2958 log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );2959 2960 break;2961 2962 case n_gappixels: // number of gap pixels2963 2964 // get the name of the output_file from the line2965 2966 sscanf(line, "%s %d", token, &ct_NGapPixels);2967 2968 log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );2969 2970 break;2971 2972 case pixel_width: // pixel width [cm]2973 2974 // get the name of the ct_file from the line2975 2976 sscanf(line, "%s %f", token, &ct_PixelWidth);2977 2978 ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));2979 2980 log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );2981 2982 break;2983 2984 case define_mirrors: // read table with the parameters of the mirrors2985 2986 log( "read_ct_file", "Table of mirrors data:\n" );2987 2988 // check whether the number of mirrors was already set2989 2990 if ( ct_NMirrors == 0 )2991 error( "read_ct_file", "NMirrors was not set.\n" );2992 2993 // allocate memory for paths list2994 2995 log( "read_ct_file", "Allocating memory for ct_data\n" );2996 2997 ct_data = new float*[ct_NMirrors];2998 2999 for (i=0; i<ct_NMirrors; i++)3000 ct_data[i] = new float[CT_NDATA];3001 3002 // read data3003 3004 log( "read_ct_file", "Reading mirrors data...\n" );3005 3006 for (i=0; i<ct_NMirrors; i++)3007 for (j=0; j<CT_NDATA; j++)3008 ctin >> ct_data[i][j];3009 3010 break;3011 3012 } // switch ( i )3013 3014 } // end while3015 3016 // end3017 3018 log( "read_ct_file", "done.\n" );3019 3020 return;3021 }3022 //!@}3023 3024 3025 //!-----------------------------------------------------------3026 2709 // @name read_QE 3027 2710 // … … 3209 2892 //!@} 3210 2893 3211 //!-----------------------------------------------------------3212 // @name read_pixels3213 //3214 // @desc read pixels data3215 //3216 // @date Fri Mar 12 16:33:34 MET 19993217 //------------------------------------------------------------3218 // @function3219 3220 //!@{3221 void3222 read_pixels(struct camera *pcam)3223 {3224 /* ifstream qefile;3225 char line[LINE_MAX_LENGTH];3226 int n, i, j, icount;3227 float qe;3228 3229 //------------------------------------------------------------3230 // first, pixels' coordinates3231 3232 pcam->inumpixels = ct_NPixels;3233 pcam->inumcentralpixels = ct_NCentralPixels;3234 pcam->inumgappixels = ct_NGapPixels;3235 pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;3236 pcam->dpixdiameter_cm = ct_PixelWidth;3237 3238 // initialize pixel numbers3239 3240 pixary = new float* [2*ct_NCentralPixels];3241 for ( i=0; i<2*ct_NCentralPixels; ++i )3242 pixary[i] = new float[2];3243 3244 pixneig = new int* [ct_NCentralPixels];3245 for ( i=0; i<ct_NCentralPixels; ++i ) {3246 pixneig[i] = new int[6];3247 for ( j=0; j<6; ++j )3248 pixneig[i][j] = -1;3249 }3250 3251 npixneig = new int[ct_NCentralPixels];3252 for ( i=0; i<ct_NCentralPixels; ++i )3253 npixneig[i] = 0;3254 3255 // generate all coordinates3256 3257 igen_pixel_coordinates(pcam);3258 3259 3260 // calculate tables of neighbours3261 3262 #ifdef __DEBUG__3263 for ( n=0 ; n<ct_NPixels ; ++n ) {3264 cout << "Para el pixel " << n << ": ";3265 for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {3266 if ( pixels_are_neig(n,i) == TRUE ) {3267 pixneig[n][npixneig[n]] = i;3268 pixneig[i][npixneig[i]] = n;3269 cout << i << ' ';3270 ++npixneig[n];3271 ++npixneig[i];3272 }3273 }3274 cout << endl << flush;3275 }3276 #else // ! __DEBUG__3277 for ( n=0 ; n<ct_NCentralPixels ; ++n )3278 for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i)3279 if ( pixels_are_neig(n,i) == TRUE ) {3280 pixneig[n][npixneig[n]] = i;3281 pixneig[i][npixneig[i]] = n;3282 ++npixneig[n];3283 ++npixneig[i];3284 }3285 #endif // ! __DEBUG__3286 3287 #ifdef __DEBUG__3288 for ( n=0 ; n<ct_NPixels ; ++n ) {3289 cout << n << ':';3290 for ( j=0; j<npixneig[n]; ++j)3291 cout << ' ' << pixneig[n][j];3292 cout << endl << flush;3293 }3294 #endif // __DEBUG__3295 3296 //------------------------------------------------------------3297 // second, pixels' QE3298 3299 // try to open the file3300 3301 log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);3302 3303 qefile.open( QE_FILE );3304 3305 // if it is wrong or does not exist, exit3306 3307 if ( qefile.bad() )3308 error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );3309 3310 // read file3311 3312 log("read_pixels", "Reading data . . .\n");3313 3314 i=-1;3315 icount = 0;3316 3317 while ( ! qefile.eof() ) {3318 3319 // get line from the file3320 3321 qefile.getline(line, LINE_MAX_LENGTH);3322 3323 // skip if comment3324 3325 if ( *line == '#' )3326 continue;3327 3328 // if it is the first valid value, it is the number of QE data points3329 3330 if ( i < 0 ) {3331 3332 // get the number of datapoints3333 3334 sscanf(line, "%d", &pointsQE);3335 3336 // allocate memory for the table of QEs3337 3338 QE = new float ** [ct_NPixels];3339 3340 for ( i=0; i<ct_NPixels; ++i ) {3341 QE[i] = new float * [2];3342 QE[i][0] = new float[pointsQE];3343 QE[i][1] = new float[pointsQE];3344 }3345 3346 QElambda = new float [pointsQE];3347 3348 for ( i=0; i<pointsQE; ++i ) {3349 qefile.getline(line, LINE_MAX_LENGTH);3350 sscanf(line, "%f", &QElambda[i]);3351 }3352 3353 i=0;3354 3355 continue;3356 }3357 3358 // get the values (num-pixel, num-datapoint, QE-value)3359 3360 if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 )3361 break;3362 3363 if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&3364 ((j-1) < pointsQE) && ((j-1) > -1) ) {3365 QE[i-1][0][j-1] = QElambda[j-1];3366 QE[i-1][1][j-1] = qe;3367 }3368 3369 if ( i > ct_NPixels) break;3370 3371 icount++;3372 3373 }3374 3375 if(icount/pointsQE < ct_NPixels){3376 error( "read_pixels", "The quantum efficiency file is faulty\n (found only %d pixels instead of %d).\n",3377 icount/pointsQE, ct_NPixels );3378 }3379 3380 // close file3381 3382 qefile.close();3383 3384 // test QE3385 3386 for(icount=0; icount< ct_NPixels; icount++){3387 for(i=0; i<pointsQE; i++){3388 if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||3389 QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){3390 error( "read_pixels", "The quantum efficiency file is faulty\n pixel %d, point %d is % f, %f\n",3391 icount, i, QE[icount][0][i], QE[icount][1][i] );3392 }3393 }3394 }3395 3396 // end3397 3398 log("read_pixels", "Done.\n");3399 */3400 }3401 //!@}3402 2894 3403 2895 //!----------------------------------------------------------- … … 3534 3026 } 3535 3027 3536 //!-----------------------------------------------------------3537 // @name pixels_are_neig3538 //3539 // @desc check whether two pixels are neighbours3540 //3541 // @var pix1 Number of the first pixel3542 // @var pix2 Number of the second pixel3543 // @return TRUE: both pixels are neighbours; FALSE: oth.3544 //3545 // @date Wed Sep 9 17:58:37 MET DST 19983546 //------------------------------------------------------------3547 // @function3548 3549 //!@{3550 int3551 pixels_are_neig(int pix1, int pix2)3552 {3553 if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +3554 SQR( pixary[pix1][1] - pixary[pix2][1] ) )3555 > ct_PixelWidth_corner_2_corner )3556 return ( FALSE );3557 else3558 return ( TRUE );3559 }3560 //!@}3561 3028 3562 3029 //!-----------------------------------------------------------
Note:
See TracChangeset
for help on using the changeset viewer.