Changeset 344 for trunk


Ignore:
Timestamp:
01/25/00 08:36:23 (25 years ago)
Author:
petry
Message:
The pixelization in previous versions was buggy.
This is the first version with a correct pixelization.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx

    r340 r344  
    2121//
    2222// $RCSfile: camera.cxx,v $
    23 // $Revision: 1.3 $
     23// $Revision: 1.4 $
    2424// $Author: petry $
    25 // $Date: 2000-01-20 18:22:17 $
     25// $Date: 2000-01-25 08:36:23 $
    2626//
    2727////////////////////////////////////////////////////////////////////////
     
    10731073          else
    10741074            inputfile.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
    1075        
     1075         
    10761076          // increase number of photons
    10771077         
    10781078          ncph++;
    1079 
     1079         
    10801080          t = cphoton.get_t() ;
    1081 
     1081         
    10821082          if(t_ini == -99999){ // this is the first photon we read from this event
    10831083            t_ini = t; // memorize time
    10841084          }
    1085 
     1085         
    10861086          // The photons don't come in chronological order!
    10871087          // Put the first photon at the center of the array by adding the constant SLICES
    1088             
     1088         
    10891089          t_chan = (int) ((t - t_ini )/ WIDTH_TIMESLICE ) + SLICES ;     
    1090 
     1090         
    10911091          if (t_chan > (2 * SLICES)){
    10921092            log(SIGNATURE, "warning, channel number (%d) exceded limit (%d).\n Setting it to %d .\n",
    1093             t_chan, (2 * SLICES), (2 * SLICES));
     1093                t_chan, (2 * SLICES), (2 * SLICES));
    10941094            t_chan = (2 * SLICES);
    10951095          }
    10961096          else if(t_chan < 0){
    10971097            log(SIGNATURE, "warning, channel number (%d) below limit (%d).\n Setting it to %d .\n",
    1098             t_chan, 0, 0);
     1098                t_chan, 0, 0);
    10991099            t_chan = 0;
    11001100          }
    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           
    11361102          // 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           
    11641104          cx = cphoton.get_x();
    11651105          cy = cphoton.get_y();
     
    11711111         
    11721112          // check if photon has valid wavelength and is inside outermost camera radius
    1173 
     1113         
    11741114          if( (wl > 800.0) || (wl < 290.0) ||
    11751115              (sqrt(cx*cx + cy*cy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)) ){
    1176            
     1116            
    11771117            // read next CPhoton
    11781118            if ( Data_From_STDIN )
     
    11831123            // go to beginning of loop, the photon is lost
    11841124            continue;
    1185 
     1125           
    11861126          }
    11871127
    11881128          // cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl;
    11891129         
    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            }
    12111137          }
    1212           else{
    1213 
    1214             nPMT = -1;
    1215 
    1216           }
    1217 
    1218           // check if outside the central camera
    1219          
    1220           if ( (nPMT < 0) || (nPMT >= ct_NCentralPixels) ) {
    1221 
    1222             // check the outer pixels
    1223             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             }
    12311138           
    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
    12691141            // read next CPhoton
    12701142            if ( Data_From_STDIN )
     
    12731145              inputfile.read ( flag, SIZE_OF_FLAGS );
    12741146           
     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           
    12751182            // go to beginning of loop
    12761183            continue;
     
    12891196         
    12901197          fnpix[nPMT] += 1.0;
    1291 
     1198         
    12921199#ifdef __ROOT__
    12931200          fnslicesum[t_chan]  += 1.0 ;
    12941201          slices[nPMT][t_chan] += 1.0 ;
    12951202#endif // __ROOT__       
    1296 
     1203         
    12971204#ifdef __DETAIL_TRIGGER__
    12981205          //
     
    13001207          //
    13011208          //
    1302 
     1209         
    13031210          Trigger.Fill( nPMT, ( t - t_ini  ) ) ;
    13041211#endif // __DETAIL_TRIGGER__
    13051212         
    13061213          // read next CPhoton
    1307 
     1214         
    13081215          if ( Data_From_STDIN )
    13091216            cin.read( flag, SIZE_OF_FLAGS );
    13101217          else
    13111218            inputfile.read ( flag, SIZE_OF_FLAGS );
    1312 
     1219         
    13131220        }  // end while, i.e. found end of event
    13141221       
    13151222        log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
    13161223            ncph, ntcph);
    1317 
     1224     
    13181225        // show number of photons
    13191226       
    13201227        //cout << ncph << " photons read . . . " << endl << flush;
    1321        
     1228     
    13221229        // skip it ?
    13231230       
     
    13941301       
    13951302        // if we should apply any kind of correction, do it here.
    1396 
     1303       
    13971304        for ( i=0; i<ct_NPixels; ++i )
    13981305          fnpix[i] *= fCorrection;
    1399 
     1306       
    14001307#ifdef __DETAIL_TRIGGER__
    1401       //   Trigger.Print() ;
    1402       cout << Trigger.Diskriminate() << endl << endl ;
     1308        //   Trigger.Print() ;
     1309        cout << Trigger.Diskriminate() << endl << endl ;
    14031310#endif // __DETAIL_TRIGGER__
    1404 
     1311       
    14051312#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           
    14441355          }
    1445        
    1446           Evt->FillPixel ( (UShort_t) i , trans ) ;
    1447          
    14481356        }
    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       
    14761383#endif // __ROOT__
    14771384       
     
    14861393        //--------------------------------------------------
    14871394       
    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        
    15221395#ifdef __TRIGGER__
    15231396       
     
    15641437        //@ If the input parameter "threshold" is 0 we find the maximum
    15651438        //@ trigger threshold this event can pass
    1566 
     1439       
    15671440        for(k=0; ( qThreshold == 0 ? (k <= iMAX_THRESHOLD_PHE) : (k == 1) ); k++){
    15681441         
     
    16661539               
    16671540                for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1; ++j ) {
    1668                  
     1541               
    16691542                  if ( fnpix[pixneig[i][j]] > q0 ) {
    16701543                   
     
    17361609                   
    17371610                    trigger = FALSE;
    1738                  
     1611                   
    17391612                  }
    17401613                 
     
    18421715         
    18431716          // calculate moments and other things
    1844 
     1717         
    18451718          moments_ptr = moments( anaPixels, fnpixclean, pixary,
    1846                               plateScale_cm2deg, 0 );
     1719                                plateScale_cm2deg, 0 );
    18471720         
    18481721          charge = moments_ptr->charge ;
     
    18611734          asymx  = moments_ptr->asymx  ;
    18621735          phiasym= moments_ptr->phi;
    1863 
     1736         
    18641737          lenwid_ptr = lenwid( anaPixels, fnpixclean, pixary,
    18651738                               plateScale_cm2deg,
    18661739                               ct_PixelWidth_corner_2_corner_half);
    1867          
    1868 
     1740         
     1741         
    18691742          // fill the diagnostic Tree
    18701743         
     
    19891862          image_data[i] =       -1.; i++;
    19901863          image_data[i] =       -1.; i++;
    1991 
     1864         
    19921865          // there should be "nvar" variables
    1993 
     1866         
    19941867          if ( i != nvar )
    19951868            error( SIGNATURE, "Wrong entry length for Ntuple.\n" );
     
    19991872         
    20001873          // put this information in the data file,
    2001 
     1874         
    20021875          if ( Write_All_Data ) {
    2003 
     1876           
    20041877            datafile << ntrigger;
    20051878            for ( i=0; i<nvar; ++i )
     
    20161889         
    20171890        } // trigger == FALSE
    2018        
     1891       
    20191892#endif // __TRIGGER__
    2020 
     1893       
    20211894        //!@' @#### Save data.
    20221895        //@'
     
    20271900        // the output file
    20281901        //--------------------------------------------------
    2029 
     1902       
    20301903        //++
    20311904        // save the image to the file
     
    20351908       
    20361909        outputfile.write( (char *)&mcevth, mcevth.mysize() );
    2037      
     1910       
    20381911#ifdef __TRIGGER__
    2039 
     1912       
    20401913        // save the image
    2041      
     1914       
    20421915        if ( (trigger == TRUE) || (Write_All_Images == TRUE) )
    20431916          outputfile.write( (char *) fnpix, ct_NPixels * sizeof( float ) );
    2044 
     1917       
    20451918#else
    2046 
     1919       
    20471920        // save the image
    2048      
     1921       
    20491922        outputfile.write( (char *) fnpix, ct_NPixels * sizeof( float ) );
    2050 
     1923       
    20511924#endif // __TRIGGER__
    2052 
     1925       
    20531926        if ( Data_From_STDIN )
    20541927          cin.read( flag, SIZE_OF_FLAGS );
     
    20571930       
    20581931      } // end while there is a next event
    2059 
     1932     
    20601933      if( !isA( flag, FLAG_END_OF_RUN   )){
    20611934        error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
     
    21011974           
    21021975          }
    2103        
     1976         
    21041977        } // end if found end of file
    21051978      } // end if found end of run
     
    21101983    } // end if else found start of run
    21111984  } // end big while loop
    2112 
    2113   //!@' @#### End of program.
    2114   //@'
    2115 
     1985 
     1986//!@' @#### End of program.
     1987//@'
     1988 
    21161989  //end my version
    21171990
    21181991#ifdef __ROOT__
    2119       //++
    2120       // put the Event to the root file
    2121       //--
    2122 
    2123       EvtTree.Write() ;
    2124       outfile.Write() ;
    2125       outfile.Close() ;
     1992  //++
     1993  // put the Event to the root file
     1994  //--
     1995 
     1996  EvtTree.Write() ;
     1997  outfile.Write() ;
     1998  outfile.Close() ;
    21261999
    21272000#endif // __ROOT__
    2128              
     2001 
    21292002  // close input file
    21302003 
     
    21342007  log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n",
    21352008       ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
    2136 
     2009 
    21372010  // close files
    21382011 
    21392012  log( SIGNATURE, "Closing files\n" );
    2140 
     2013 
    21412014  inputfile.close();
    21422015  outputfile.close();
     
    27972670    /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
    27982671
    2799     pcam->dpixsizefactor[ipixno] = 1.;
     2672    pcam->dpixsizefactor[ipixno-1] = 1.;
    28002673
    28012674    in = 0;
     
    31042977//
    31052978// $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//
    31062986// Revision 1.2  1999/11/19 08:40:42  harald
    31072987// Now it is possible to compile the camera programm under osf1.
Note: See TracChangeset for help on using the changeset viewer.