Changeset 1512 for trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
- Timestamp:
- 09/04/02 10:57:42 (22 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
r1419 r1512 21 21 // 22 22 // $RCSfile: camera.cxx,v $ 23 // $Revision: 1.4 0$23 // $Revision: 1.41 $ 24 24 // $Author: blanch $ 25 // $Date: 2002-0 7-16 16:15:22 $25 // $Date: 2002-09-04 09:57:42 $ 26 26 // 27 27 //////////////////////////////////////////////////////////////////////// … … 72 72 #include "MMcTrigHeader.hxx" 73 73 #include "MMcFadcHeader.hxx" 74 #include "MGeomCamMagic.h" 75 #include "MGeomPix.h" 74 76 75 77 /*!@" … … 495 497 UInt_t corsika = 5200 ; 496 498 497 498 struct camera cam; // structure holding the camera definition 499 MGeomCamMagic camgeom; // structure holding the camera definition 499 500 500 501 … … 759 760 760 761 read_ct_file(); 761 762 // read camera setup763 764 read_pixels(&cam);765 762 766 763 Int_t Lev0, Lev1; … … 1112 1109 1113 1110 k = produce_nsbrates( starfieldname, 1114 &cam ,1111 &camgeom, 1115 1112 nsbrate_phepns); 1116 1113 … … 1122 1119 // calculate diffuse rate correcting for the pixel size 1123 1120 1124 for(i=0; i<cam.inumpixels; i++){ 1125 diffnsb_phepns[i] = meanNSB * 1126 cam.dpixsizefactor[i] * cam.dpixsizefactor[i]; 1121 const double size_inner = camgeom[0].GetR(); 1122 1123 for(UInt_t ui=0; ui<camgeom.GetNumPixels(); ui++){ 1124 const double actual_size = camgeom[ui].GetR(); 1125 diffnsb_phepns[ui] = meanNSB * 1126 actual_size*actual_size/(size_inner*size_inner); 1127 1127 } 1128 1128 … … 1134 1134 zenfactor = pow(10., -0.4 * ext[j] ); 1135 1135 1136 for( i=0; i<cam.inumpixels;i++){1137 nsb_phepns[ i]+=diffnsb_phepns[i]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[i][j];1136 for(UInt_t ui=0; ui<camgeom.GetNumPixels();ui++){ 1137 nsb_phepns[ui]+=diffnsb_phepns[ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j]; 1138 1138 } 1139 1139 } … … 1327 1327 1328 1328 k = produce_phes( inputfile, 1329 &cam ,1329 &camgeom, 1330 1330 WAVEBANDBOUND1, 1331 1331 WAVEBANDBOUND6, … … 1349 1349 1350 1350 // Fill trigger and fadc response in the trigger class from the database 1351 for( i=0;i<cam.inumpixels;i++)1352 if(nsb_phepns[ i]>0.0){1353 k=lons.GetResponse(nsb_phepns[ i],0.1,1351 for(UInt_t ui=0;ui<camgeom.GetNumPixels();ui++) 1352 if(nsb_phepns[ui]>0.0){ 1353 k=lons.GetResponse(nsb_phepns[ui],0.1, 1354 1354 & nsb_trigresp[0],& nsb_fadcresp[0]); 1355 1355 if(k==0){ … … 1357 1357 exit(1); 1358 1358 } 1359 Trigger.AddNSB( i,nsb_trigresp);1360 fadc.AddSignal( i,nsb_fadcresp);1359 Trigger.AddNSB(ui,nsb_trigresp); 1360 fadc.AddSignal(ui,nsb_fadcresp); 1361 1361 } 1362 1362 }// end if(simulateNSB && nphe2NSB<=inumphe) ... … … 1387 1387 cout << "Total number of phes: " << inumphe<<" + "; 1388 1388 inumphensb=0; 1389 for ( i=0;i<cam.inumpixels;i++){1390 inumphensb+=nsb_phepns[ i]*TOTAL_TRIGGER_TIME;1389 for (UInt_t ui=0;ui<camgeom.GetNumPixels();ui++){ 1390 inumphensb+=nsb_phepns[ui]*TOTAL_TRIGGER_TIME; 1391 1391 } 1392 1392 cout<<inumphensb<<endl; … … 2751 2751 2752 2752 /******** bpoint_is_in_pix() ***************************************/ 2753 /*2754 int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){2755 // return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system)2756 // the pixel is assumed to be a "closed set"2757 2758 double a, b; // a = length of one of the edges of one pixel, b = half the width of one pixel2759 double c, xx, yy; // auxiliary variable2760 2761 b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];2762 a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];2763 c = 1./sqrt(3.);2764 2765 if((ipixnum < 0)||(ipixnum >= pcam->inumpixels)){2766 fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);2767 fprintf(stderr, "Exiting.\n");2768 exit(203);2769 }2770 xx = dx - pcam->dxc[ipixnum];2771 yy = dy - pcam->dyc[ipixnum];2772 2773 if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) ||2774 ((0. < xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a))) ){2775 return(TRUE); // point is inside2776 }2777 else{2778 return(FALSE); // point is outside2779 }2780 }2781 */2782 2753 2783 2754 #define sqrt13 0.577350269 // = 1./sqrt(3.) 2784 2755 #define sqrt3 1.732050807 // = sqrt(3.) 2785 2756 2786 int bpoint_is_in_pix(double dx, double dy, struct camera *pcam)2757 int bpoint_is_in_pix(double dx, double dy, MGeomCamMagic *pgeomcam) 2787 2758 { 2788 2759 /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */ … … 2794 2765 */ 2795 2766 2796 const int num = pcam->inumpixels; 2797 const double diam = pcam->dpixdiameter_cm; 2798 2799 const double diam2 = diam/2; 2800 const double diam3 = diam/sqrt3; 2801 2802 for (int i=0; i<num; i++) 2767 const int numN = pgeomcam->GetNumPixels(); 2768 2769 for (int i=0; i<numN; i++) 2803 2770 { 2804 const double size = pcam->dpixsizefactor[i]; 2805 2806 const double b = diam2 * size; 2807 const double a = diam3 * size; 2808 2809 const double xx = dx - pcam->dxc[i]; 2810 const double yy = dy - pcam->dyc[i]; 2811 2812 if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) || 2813 ((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){ 2814 2815 return i; // inside i 2771 MGeomPix &pixel = (*pgeomcam)[i]; 2772 const double b = pixel.GetR()/2; 2773 const double a = pixel.GetR()/sqrt3; 2774 2775 const double xx = dx - pixel.GetX(); 2776 const double yy = dy - pixel.GetY(); 2777 2778 if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) || 2779 ((0. < xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a))) ){ 2780 2781 return i; // inside i 2816 2782 } 2817 2783 2818 2784 // return -1; // outside 2819 2785 } 2786 2820 2787 return -1; // outside 2821 2788 } … … 2917 2884 2918 2885 int produce_phes( FILE *sp, // the input file 2919 struct camera *cam, // the camera layout2886 MGeomCamMagic *camgeom, // the camera layout 2920 2887 float minwl_nm, // the minimum accepted wavelength 2921 2888 float maxwl_nm, // the maximum accepted wavelength … … 2930 2897 ){ 2931 2898 2932 static int i, k, ipixnum; 2899 static uint i; 2900 static int k, ipixnum; 2933 2901 static float cx, cy, wl, qe, t; 2934 2902 static MCCphoton photon; 2935 2903 static float **qept; 2936 2904 static char flag[SIZE_OF_FLAGS + 1]; 2937 static float radius ;2905 static float radius_mm; 2938 2906 static float SFR , C3 , C2 , C1; 2939 2907 const double pi = 3.14159; … … 2941 2909 // reset variables 2942 2910 2943 for ( i=0; i<cam ->inumpixels; ++i ){2911 for ( i=0; i<camgeom->GetNumPixels(); ++i ){ 2944 2912 2945 2913 nphe[i] = 0.0; … … 2949 2917 *incph = 0; 2950 2918 2951 radius = cam->dxc[cam->inumpixels-1]2952 + 1.5*cam->dpixdiameter_cm*cam->dpixsizefactor[cam->inumpixels-1]; 2919 radius_mm = camgeom->GetMaxRadius(); 2920 2953 2921 2954 2922 //- - - - - - - - - - - - - - - - - - - - - - - - - … … 3043 3011 // check if photon has valid wavelength and is inside outermost camera radius 3044 3012 3045 if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius) ){3013 if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy)*10 > radius_mm ) ){ 3046 3014 continue; 3047 3015 3048 3016 } 3049 3017 3050 ipixnum = bpoint_is_in_pix(cx , cy, cam);3018 ipixnum = bpoint_is_in_pix(cx*10, cy*10, camgeom); 3051 3019 3052 3020 // -1 = the photon is in none of the pixels … … 3121 3089 3122 3090 int produce_nsbrates( char *iname, // the starfield input file name 3123 struct camera *cam, // camera layout3091 MGeomCamMagic *camgeom, // camera layout 3124 3092 float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function: 3125 3093 // the NSB rates in phe/ns for each pixel 3126 3094 ){ 3127 3095 3128 int i, j, k, ii; // counters 3096 uint i, j; // counters 3097 int k, ii; // counters 3129 3098 3130 3099 MTrigger trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm); … … 3212 3181 // initialize the rate array 3213 3182 3214 for(j = 0; j<cam ->inumpixels; j++){ // loop over pixels3183 for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels 3215 3184 rate_phepns[j][i] = 0.; 3216 3185 } … … 3231 3200 3232 3201 k = produce_phes( infile, 3233 cam ,3202 camgeom, 3234 3203 wl_nm[i], 3235 3204 wl_nm[i+1], … … 3248 3217 } 3249 3218 3250 for(j = 0; j<cam ->inumpixels; j++){ // loop over pixels3219 for(j = 0; j<camgeom->GetNumPixels(); j++){ // loop over pixels 3251 3220 rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS; 3252 3221 } … … 3409 3378 // 3410 3379 // $Log: not supported by cvs2svn $ 3380 // Revision 1.40 2002/07/16 16:15:22 blanch 3381 // A first implementation of the Star field rotation has been done. 3382 // 3411 3383 // Revision 1.39 2002/04/27 10:48:39 blanch 3412 3384 // Some unused varibles have been removed. … … 3528 3500 // 3529 3501 // $Log: not supported by cvs2svn $ 3502 // Revision 1.40 2002/07/16 16:15:22 blanch 3503 // A first implementation of the Star field rotation has been done. 3504 // 3530 3505 // Revision 1.39 2002/04/27 10:48:39 blanch 3531 3506 // Some unused varibles have been removed.
Note:
See TracChangeset
for help on using the changeset viewer.