Changeset 344 for trunk/MagicSoft/Simulation/Detector/Camera
- Timestamp:
- 01/25/00 08:36:23 (25 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
r340 r344 21 21 // 22 22 // $RCSfile: camera.cxx,v $ 23 // $Revision: 1. 3$23 // $Revision: 1.4 $ 24 24 // $Author: petry $ 25 // $Date: 2000-01-2 0 18:22:17$25 // $Date: 2000-01-25 08:36:23 $ 26 26 // 27 27 //////////////////////////////////////////////////////////////////////// … … 1073 1073 else 1074 1074 inputfile.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS ); 1075 1075 1076 1076 // increase number of photons 1077 1077 1078 1078 ncph++; 1079 1079 1080 1080 t = cphoton.get_t() ; 1081 1081 1082 1082 if(t_ini == -99999){ // this is the first photon we read from this event 1083 1083 t_ini = t; // memorize time 1084 1084 } 1085 1085 1086 1086 // The photons don't come in chronological order! 1087 1087 // Put the first photon at the center of the array by adding the constant SLICES 1088 1088 1089 1089 t_chan = (int) ((t - t_ini )/ WIDTH_TIMESLICE ) + SLICES ; 1090 1090 1091 1091 if (t_chan > (2 * SLICES)){ 1092 1092 log(SIGNATURE, "warning, channel number (%d) exceded limit (%d).\n Setting it to %d .\n", 1093 1093 t_chan, (2 * SLICES), (2 * SLICES)); 1094 1094 t_chan = (2 * SLICES); 1095 1095 } 1096 1096 else if(t_chan < 0){ 1097 1097 log(SIGNATURE, "warning, channel number (%d) below limit (%d).\n Setting it to %d .\n", 1098 1098 t_chan, 0, 0); 1099 1099 t_chan = 0; 1100 1100 } 1101 /*!@' 1102 1103 @#### Pixelization (for the central pixels). 1104 1105 In order to calculate the coordinates, we use the 1106 change of system described in the documentation 1107 of the source code of |pixel\_coord.cxx|. 1108 Then, we will use simply the matrix of change 1109 from one system to the other. In our case, this is: 1110 1111 @[ 1112 \begin{bmatrix}X\\Y\\\end{bmatrix} 1113 = 1114 \begin{bmatrix} 1115 1 & \cos(60^\circ)\\ 1116 0 & \sin(60^\circ)\\ 1117 \end{bmatrix} 1118 \begin{bmatrix}I\\J\\\end{bmatrix} 1119 @] 1120 1121 and hence 1122 1123 @[ 1124 \begin{bmatrix}I\\J\\\end{bmatrix} 1125 = 1126 \begin{bmatrix} 1127 1 & -\frac{\cos(60^\circ)}{\sin(60^\circ)}\\ 1128 0 &\frac{1}{\sin(60^\circ)}\\ 1129 \end{bmatrix} 1130 \begin{bmatrix}X\\Y\\\end{bmatrix} 1131 @] 1132 1133 */ 1134 1135 //+++ 1101 1136 1102 // Pixelization 1137 //--- 1138 1139 // calculate ij-coordinates 1140 1141 // We use a change of coordinate system, using the following 1142 // matrix of change (m^-1) (this is taken from Mathematica output). 1143 /* 1144 * In[1]:= m={{1,cos60},{0,sin60}}; MatrixForm[m] 1145 * 1146 * Out[1]//MatrixForm= 1 cos60 1147 * 1148 * 0 sin60 1149 * 1150 * In[2]:= inv=Inverse[m]; MatrixForm[inv] 1151 * 1152 * Out[2]//MatrixForm= cos60 1153 * -(-----) 1154 * 1 sin60 1155 * 1156 * 1 1157 * ----- 1158 * 0 sin60 1159 * 1160 */ 1161 1162 // go to IJ-coordinate system 1163 1103 1164 1104 cx = cphoton.get_x(); 1165 1105 cy = cphoton.get_y(); … … 1171 1111 1172 1112 // check if photon has valid wavelength and is inside outermost camera radius 1173 1113 1174 1114 if( (wl > 800.0) || (wl < 290.0) || 1175 1115 (sqrt(cx*cx + cy*cy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)) ){ 1176 1116 1177 1117 // read next CPhoton 1178 1118 if ( Data_From_STDIN ) … … 1183 1123 // go to beginning of loop, the photon is lost 1184 1124 continue; 1185 1125 1186 1126 } 1187 1127 1188 1128 // cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl; 1189 1129 1190 ci = floor( (cx - cy*COS60/SIN60)/ ct_2Apot + 0.5); 1191 cj = floor( (cy/SIN60) / ct_2Apot + 0.5); 1192 1193 ici = (int)(ci); 1194 icj = (int)(cj); 1195 1196 iici = ici+PIX_ARRAY_HALF_SIDE; 1197 iicj = icj+PIX_ARRAY_HALF_SIDE; 1198 1199 // is it inside the array? 1200 1201 if ( (iici > 0) && (iici < PIX_ARRAY_SIDE) && 1202 (iicj > 0) && (iicj < PIX_ARRAY_SIDE) ) { 1203 1204 // try to put into pixel 1205 1206 // obtain the pixel number for this photon 1207 1208 nPMT = (int) 1209 pixels[ici+PIX_ARRAY_HALF_SIDE][icj+PIX_ARRAY_HALF_SIDE][PIXNUM]; 1210 1130 nPMT = -1; 1131 1132 for(i=0; i<ct_NPixels; i++){ 1133 if( bpoint_is_in_pix( cx, cy, i, &cam) ){ 1134 nPMT = i; 1135 break; 1136 } 1211 1137 } 1212 else{1213 1214 nPMT = -1;1215 1216 }1217 1218 // check if outside the central camera1219 1220 if ( (nPMT < 0) || (nPMT >= ct_NCentralPixels) ) {1221 1222 // check the outer pixels1223 nPMT = -1;1224 1225 for(i=ct_NCentralPixels; i<ct_NPixels; i++){1226 if( bpoint_is_in_pix( cx, cy, i, &cam) ){1227 nPMT = i;1228 break;1229 }1230 }1231 1138 1232 if(nPMT==-1){// the photon is in none of the pixels 1233 1234 // read next CPhoton 1235 if ( Data_From_STDIN ) 1236 cin.read( flag, SIZE_OF_FLAGS ); 1237 else 1238 inputfile.read ( flag, SIZE_OF_FLAGS ); 1239 1240 // go to beginning of loop, the photon is lost 1241 continue; 1242 } 1243 1244 } 1245 1246 #ifdef __QE__ 1247 1248 //!@' @#### QE simulation. 1249 //@' 1250 1251 //+++ 1252 // QE simulation 1253 //--- 1254 1255 // find data point to be used in Lagrange interpolation (-> k) 1256 1257 qeptr = (float **)QE[nPMT]; 1258 1259 FindLagrange(qeptr,k,wl); 1260 1261 // if random > quantum efficiency, reject it 1262 1263 qe = Lagrange(qeptr,k,wl) / 100.0; 1264 1265 // fprintf(stdout, "%f\n", qe); 1266 1267 if ( RandomNumber > qe ) { 1268 1139 if(nPMT==-1){// the photon is in none of the pixels 1140 1269 1141 // read next CPhoton 1270 1142 if ( Data_From_STDIN ) … … 1273 1145 inputfile.read ( flag, SIZE_OF_FLAGS ); 1274 1146 1147 // go to beginning of loop, the photon is lost 1148 continue; 1149 } 1150 1151 1152 1153 #ifdef __QE__ 1154 1155 //!@' @#### QE simulation. 1156 //@' 1157 1158 //+++ 1159 // QE simulation 1160 //--- 1161 1162 // find data point to be used in Lagrange interpolation (-> k) 1163 1164 qeptr = (float **)QE[nPMT]; 1165 1166 FindLagrange(qeptr,k,wl); 1167 1168 // if random > quantum efficiency, reject it 1169 1170 qe = Lagrange(qeptr,k,wl) / 100.0; 1171 1172 // fprintf(stdout, "%f\n", qe); 1173 1174 if ( RandomNumber > qe ) { 1175 1176 // read next CPhoton 1177 if ( Data_From_STDIN ) 1178 cin.read( flag, SIZE_OF_FLAGS ); 1179 else 1180 inputfile.read ( flag, SIZE_OF_FLAGS ); 1181 1275 1182 // go to beginning of loop 1276 1183 continue; … … 1289 1196 1290 1197 fnpix[nPMT] += 1.0; 1291 1198 1292 1199 #ifdef __ROOT__ 1293 1200 fnslicesum[t_chan] += 1.0 ; 1294 1201 slices[nPMT][t_chan] += 1.0 ; 1295 1202 #endif // __ROOT__ 1296 1203 1297 1204 #ifdef __DETAIL_TRIGGER__ 1298 1205 // … … 1300 1207 // 1301 1208 // 1302 1209 1303 1210 Trigger.Fill( nPMT, ( t - t_ini ) ) ; 1304 1211 #endif // __DETAIL_TRIGGER__ 1305 1212 1306 1213 // read next CPhoton 1307 1214 1308 1215 if ( Data_From_STDIN ) 1309 1216 cin.read( flag, SIZE_OF_FLAGS ); 1310 1217 else 1311 1218 inputfile.read ( flag, SIZE_OF_FLAGS ); 1312 1219 1313 1220 } // end while, i.e. found end of event 1314 1221 1315 1222 log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n", 1316 1223 ncph, ntcph); 1317 1224 1318 1225 // show number of photons 1319 1226 1320 1227 //cout << ncph << " photons read . . . " << endl << flush; 1321 1228 1322 1229 // skip it ? 1323 1230 … … 1394 1301 1395 1302 // if we should apply any kind of correction, do it here. 1396 1303 1397 1304 for ( i=0; i<ct_NPixels; ++i ) 1398 1305 fnpix[i] *= fCorrection; 1399 1306 1400 1307 #ifdef __DETAIL_TRIGGER__ 1401 1402 1308 // Trigger.Print() ; 1309 cout << Trigger.Diskriminate() << endl << endl ; 1403 1310 #endif // __DETAIL_TRIGGER__ 1404 1311 1405 1312 #ifdef __ROOT__ 1406 1407 // 1408 // Fill the header of this event 1409 // 1410 1411 Evt->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ; 1412 1413 // now put out the data of interest 1414 // 1415 // 1. -> look for the first slice with signal 1416 // 1417 1418 for ( i=0; i<(2 * SLICES); i++ ) 1419 if ( fnslicesum[i] > 0. ) 1420 break ; 1421 1422 startchan = i ; 1423 1424 // 1425 // copy the slices out of the big array 1426 // 1427 // put the first slice with signal to slice 4 1428 // 1429 1430 for (i=0; i<ct_NPixels; i++ ) 1431 for ( ii=(startchan-3); ii < (startchan+12); ii++ ) 1432 slices2 [i][ii-startchan+3] = slices [i][ii] ; 1433 1434 1435 // 1436 // if a pixes has a signal put it to the MRawEvt 1437 // 1438 1439 for (i=0; i<ct_NPixels; i++ ) { 1440 if ( fnpix[i] > 0 ) { 1441 1442 for ( ii=0; ii < 15; ii++ ) { 1443 trans [ii] = slices2[i][ii] ; 1313 1314 // 1315 // Fill the header of this event 1316 // 1317 1318 Evt->FillHeader ( (UShort_t) (ntshow + nshow) , 20 ) ; 1319 1320 // now put out the data of interest 1321 // 1322 // 1. -> look for the first slice with signal 1323 // 1324 1325 for ( i=0; i<(2 * SLICES); i++ ) 1326 if ( fnslicesum[i] > 0. ) 1327 break ; 1328 1329 startchan = i ; 1330 1331 // 1332 // copy the slices out of the big array 1333 // 1334 // put the first slice with signal to slice 4 1335 // 1336 1337 for (i=0; i<ct_NPixels; i++ ) 1338 for ( ii=(startchan-3); ii < (startchan+12); ii++ ) 1339 slices2 [i][ii-startchan+3] = slices [i][ii] ; 1340 1341 1342 // 1343 // if a pixes has a signal put it to the MRawEvt 1344 // 1345 1346 for (i=0; i<ct_NPixels; i++ ) { 1347 if ( fnpix[i] > 0 ) { 1348 1349 for ( ii=0; ii < 15; ii++ ) { 1350 trans [ii] = slices2[i][ii] ; 1351 } 1352 1353 Evt->FillPixel ( (UShort_t) i , trans ) ; 1354 1444 1355 } 1445 1446 Evt->FillPixel ( (UShort_t) i , trans ) ;1447 1448 1356 } 1449 } 1450 1451 // 1452 // 1453 // 1454 1455 McEvt->Fill( (UShort_t) mcevth.get_primary() , 1456 mcevth.get_energy(), 1457 mcevth.get_theta(), 1458 mcevth.get_phi(), 1459 mcevth.get_core(), 1460 mcevth.get_coreX(), 1461 mcevth.get_coreY(), 1462 flli, 1463 ulli, ulli, ulli, ulli, ulli ) ; 1464 1465 // 1466 // write it out to the file outfile 1467 // 1468 1469 EvtTree.Fill() ; 1470 1471 // clear all 1472 1473 Evt->Clear() ; 1474 McEvt->Clear() ; 1475 1357 1358 // 1359 // 1360 // 1361 1362 McEvt->Fill( (UShort_t) mcevth.get_primary() , 1363 mcevth.get_energy(), 1364 mcevth.get_theta(), 1365 mcevth.get_phi(), 1366 mcevth.get_core(), 1367 mcevth.get_coreX(), 1368 mcevth.get_coreY(), 1369 flli, 1370 ulli, ulli, ulli, ulli, ulli ) ; 1371 1372 // 1373 // write it out to the file outfile 1374 // 1375 1376 EvtTree.Fill() ; 1377 1378 // clear all 1379 1380 Evt->Clear() ; 1381 McEvt->Clear() ; 1382 1476 1383 #endif // __ROOT__ 1477 1384 … … 1486 1393 //-------------------------------------------------- 1487 1394 1488 #ifdef __DEBUG__1489 printf("\n");1490 1491 for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {1492 1493 for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {1494 1495 if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {1496 1497 if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {1498 1499 printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow,1500 (int)pixels[ici][icj][PIXNUM],1501 pixels[ici][icj][PIXX],1502 pixels[ici][icj][PIXY],1503 fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);1504 1505 }1506 1507 }1508 1509 }1510 1511 }1512 1513 for (i=0; i<ct_NPixels; ++i) {1514 printf("%d (%d): ", i, npixneig[i]);1515 for (j=0; j<npixneig[i]; ++i)1516 printf(" %d", pixneig[i][j]);1517 printf("\n");1518 }1519 1520 #endif // __DEBUG__1521 1522 1395 #ifdef __TRIGGER__ 1523 1396 … … 1564 1437 //@ If the input parameter "threshold" is 0 we find the maximum 1565 1438 //@ trigger threshold this event can pass 1566 1439 1567 1440 for(k=0; ( qThreshold == 0 ? (k <= iMAX_THRESHOLD_PHE) : (k == 1) ); k++){ 1568 1441 … … 1666 1539 1667 1540 for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1; ++j ) { 1668 1541 1669 1542 if ( fnpix[pixneig[i][j]] > q0 ) { 1670 1543 … … 1736 1609 1737 1610 trigger = FALSE; 1738 1611 1739 1612 } 1740 1613 … … 1842 1715 1843 1716 // calculate moments and other things 1844 1717 1845 1718 moments_ptr = moments( anaPixels, fnpixclean, pixary, 1846 1719 plateScale_cm2deg, 0 ); 1847 1720 1848 1721 charge = moments_ptr->charge ; … … 1861 1734 asymx = moments_ptr->asymx ; 1862 1735 phiasym= moments_ptr->phi; 1863 1736 1864 1737 lenwid_ptr = lenwid( anaPixels, fnpixclean, pixary, 1865 1738 plateScale_cm2deg, 1866 1739 ct_PixelWidth_corner_2_corner_half); 1867 1868 1740 1741 1869 1742 // fill the diagnostic Tree 1870 1743 … … 1989 1862 image_data[i] = -1.; i++; 1990 1863 image_data[i] = -1.; i++; 1991 1864 1992 1865 // there should be "nvar" variables 1993 1866 1994 1867 if ( i != nvar ) 1995 1868 error( SIGNATURE, "Wrong entry length for Ntuple.\n" ); … … 1999 1872 2000 1873 // put this information in the data file, 2001 1874 2002 1875 if ( Write_All_Data ) { 2003 1876 2004 1877 datafile << ntrigger; 2005 1878 for ( i=0; i<nvar; ++i ) … … 2016 1889 2017 1890 } // trigger == FALSE 2018 1891 2019 1892 #endif // __TRIGGER__ 2020 1893 2021 1894 //!@' @#### Save data. 2022 1895 //@' … … 2027 1900 // the output file 2028 1901 //-------------------------------------------------- 2029 1902 2030 1903 //++ 2031 1904 // save the image to the file … … 2035 1908 2036 1909 outputfile.write( (char *)&mcevth, mcevth.mysize() ); 2037 1910 2038 1911 #ifdef __TRIGGER__ 2039 1912 2040 1913 // save the image 2041 1914 2042 1915 if ( (trigger == TRUE) || (Write_All_Images == TRUE) ) 2043 1916 outputfile.write( (char *) fnpix, ct_NPixels * sizeof( float ) ); 2044 1917 2045 1918 #else 2046 1919 2047 1920 // save the image 2048 1921 2049 1922 outputfile.write( (char *) fnpix, ct_NPixels * sizeof( float ) ); 2050 1923 2051 1924 #endif // __TRIGGER__ 2052 1925 2053 1926 if ( Data_From_STDIN ) 2054 1927 cin.read( flag, SIZE_OF_FLAGS ); … … 2057 1930 2058 1931 } // end while there is a next event 2059 1932 2060 1933 if( !isA( flag, FLAG_END_OF_RUN )){ 2061 1934 error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag ); … … 2101 1974 2102 1975 } 2103 1976 2104 1977 } // end if found end of file 2105 1978 } // end if found end of run … … 2110 1983 } // end if else found start of run 2111 1984 } // end big while loop 2112 2113 2114 2115 1985 1986 //!@' @#### End of program. 1987 //@' 1988 2116 1989 //end my version 2117 1990 2118 1991 #ifdef __ROOT__ 2119 2120 2121 2122 2123 2124 2125 1992 //++ 1993 // put the Event to the root file 1994 //-- 1995 1996 EvtTree.Write() ; 1997 outfile.Write() ; 1998 outfile.Close() ; 2126 1999 2127 2000 #endif // __ROOT__ 2128 2001 2129 2002 // close input file 2130 2003 … … 2134 2007 log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n", 2135 2008 ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow); 2136 2009 2137 2010 // close files 2138 2011 2139 2012 log( SIGNATURE, "Closing files\n" ); 2140 2013 2141 2014 inputfile.close(); 2142 2015 outputfile.close(); … … 2797 2670 /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */ 2798 2671 2799 pcam->dpixsizefactor[ipixno ] = 1.;2672 pcam->dpixsizefactor[ipixno-1] = 1.; 2800 2673 2801 2674 in = 0; … … 3104 2977 // 3105 2978 // $Log: not supported by cvs2svn $ 2979 // Revision 1.3 2000/01/20 18:22:17 petry 2980 // Found little bug which makes camera crash if it finds a photon 2981 // of invalid wavelength. This bug is now fixed and the range 2982 // of valid wavelengths extended to 290 - 800 nm. 2983 // This is in preparation for the NSB simulation to come. 2984 // Dirk 2985 // 3106 2986 // Revision 1.2 1999/11/19 08:40:42 harald 3107 2987 // Now it is possible to compile the camera programm under osf1.
Note:
See TracChangeset
for help on using the changeset viewer.
