Changeset 344 for trunk/MagicSoft/Simulation

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:

• trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (modified) (25 diffs)

trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx

@@ -21 +21 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.3 $
+// $Revision: 1.4 $
 // $Author: petry $ 
-// $Date: 2000-01-20 18:22:17 $
+// $Date: 2000-01-25 08:36:23 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -1073 +1073 @@
           else
             inputfile.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
-        
+          
           // increase number of photons
           
           ncph++;
-
+          
           t = cphoton.get_t() ; 
-
+          
           if(t_ini == -99999){ // this is the first photon we read from this event
             t_ini = t; // memorize time
           }
-
+          
           // The photons don't come in chronological order!
           // Put the first photon at the center of the array by adding the constant SLICES
-             
+          
           t_chan = (int) ((t - t_ini )/ WIDTH_TIMESLICE ) + SLICES ;      
-
+          
           if (t_chan > (2 * SLICES)){
             log(SIGNATURE, "warning, channel number (%d) exceded limit (%d).\n Setting it to %d .\n",
-            t_chan, (2 * SLICES), (2 * SLICES));
+                t_chan, (2 * SLICES), (2 * SLICES));
             t_chan = (2 * SLICES);
           }
           else if(t_chan < 0){
             log(SIGNATURE, "warning, channel number (%d) below limit (%d).\n Setting it to %d .\n",
-            t_chan, 0, 0);
+                t_chan, 0, 0);
             t_chan = 0;
           } 
-          /*!@'
-            
-            @#### Pixelization (for the central pixels).
-            
-            In order to calculate the coordinates, we use the
-            change of system described in the documentation
-            of the source code of |pixel\_coord.cxx|.
-            Then, we will use simply the matrix of change
-            from one system to the other. In our case, this is:
-            
-            @[
-            \begin{bmatrix}X\\Y\\\end{bmatrix}                                
-            =
-            \begin{bmatrix}
-            1 & \cos(60^\circ)\\
-            0 & \sin(60^\circ)\\
-            \end{bmatrix}                                
-            \begin{bmatrix}I\\J\\\end{bmatrix}                                
-            @]
-            
-            and hence
-            
-            @[
-            \begin{bmatrix}I\\J\\\end{bmatrix}                                
-            =
-            \begin{bmatrix}    
-            1 & -\frac{\cos(60^\circ)}{\sin(60^\circ)}\\
-            0 &\frac{1}{\sin(60^\circ)}\\
-            \end{bmatrix}                                
-            \begin{bmatrix}X\\Y\\\end{bmatrix}                                
-            @]
-            
-          */
-          
-          //+++
+            
           // Pixelization
-          //---
-          
-          // calculate ij-coordinates
-          
-          // We use a change of coordinate system, using the following 
-          // matrix of change (m^-1) (this is taken from Mathematica output).
-          /*
-           * In[1]:= m={{1,cos60},{0,sin60}}; MatrixForm[m]
-           *
-           * Out[1]//MatrixForm= 1       cos60
-           * 
-           *                     0       sin60
-           * 
-           * In[2]:= inv=Inverse[m]; MatrixForm[inv]
-           * 
-           * Out[2]//MatrixForm=              cos60
-           *                                -(-----)
-           *                       1          sin60
-           * 
-           *                                    1
-           *                                  -----
-           *                       0          sin60
-           * 
-           */
-          
-          // go to IJ-coordinate system
-          
+            
           cx = cphoton.get_x();
           cy = cphoton.get_y(); 
@@ -1171 +1111 @@
           
           // check if photon has valid wavelength and is inside outermost camera radius
-
+          
           if( (wl > 800.0) || (wl < 290.0) ||
               (sqrt(cx*cx + cy*cy) > (cam.dxc[ct_NPixels-1]+1.5*ct_PixelWidth)) ){ 
-           
+            
             // read next CPhoton
             if ( Data_From_STDIN ) 
@@ -1183 +1123 @@
             // go to beginning of loop, the photon is lost
             continue;
-
+            
           }
 
           // cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl;
           
-          ci = floor( (cx - cy*COS60/SIN60)/ ct_2Apot + 0.5);
-          cj = floor( (cy/SIN60) / ct_2Apot + 0.5);
-          
-          ici = (int)(ci);
-          icj = (int)(cj);
-          
-          iici = ici+PIX_ARRAY_HALF_SIDE;
-          iicj = icj+PIX_ARRAY_HALF_SIDE;
-          
-          // is it inside the array?
-          
-          if ( (iici > 0) && (iici < PIX_ARRAY_SIDE) &&
-               (iicj > 0) && (iicj < PIX_ARRAY_SIDE) ) {
-                  
-            // try to put into pixel
-            
-            // obtain the pixel number for this photon
-            
-            nPMT = (int)
-              pixels[ici+PIX_ARRAY_HALF_SIDE][icj+PIX_ARRAY_HALF_SIDE][PIXNUM];
-          
+          nPMT = -1;
+
+          for(i=0; i<ct_NPixels; i++){
+            if( bpoint_is_in_pix( cx, cy, i, &cam) ){
+              nPMT = i;
+              break;
+            }
           }
-          else{
-
-            nPMT = -1;
-
-          }
-
-          // check if outside the central camera
-          
-          if ( (nPMT < 0) || (nPMT >= ct_NCentralPixels) ) {
-
-            // check the outer pixels
-            nPMT = -1;
-
-            for(i=ct_NCentralPixels; i<ct_NPixels; i++){
-              if( bpoint_is_in_pix( cx, cy, i, &cam) ){
-                nPMT = i;
-                break;
-              }
-            }
            
-            if(nPMT==-1){// the photon is in none of the pixels
-
-              // read next CPhoton
-              if ( Data_From_STDIN ) 
-                cin.read( flag, SIZE_OF_FLAGS );
-              else
-                inputfile.read ( flag, SIZE_OF_FLAGS );
-              
-              // go to beginning of loop, the photon is lost
-              continue;
-            }
-            
-          }
-          
-#ifdef __QE__
-          
-          //!@' @#### QE simulation.
-          //@'
-          
-          //+++
-          // QE simulation
-          //---
-          
-          // find data point to be used in Lagrange interpolation (-> k)
-          
-          qeptr = (float **)QE[nPMT];
-          
-          FindLagrange(qeptr,k,wl);
-          
-          // if random > quantum efficiency, reject it
-          
-          qe = Lagrange(qeptr,k,wl) / 100.0;
-
-          // fprintf(stdout, "%f\n", qe);
-          
-          if ( RandomNumber > qe ) {
-            
+          if(nPMT==-1){// the photon is in none of the pixels
+
             // read next CPhoton
             if ( Data_From_STDIN ) 
@@ -1273 +1145 @@
               inputfile.read ( flag, SIZE_OF_FLAGS );
             
+            // go to beginning of loop, the photon is lost
+            continue;
+          }
+          
+        
+          
+#ifdef __QE__
+          
+          //!@' @#### QE simulation.
+          //@'
+          
+          //+++
+          // QE simulation
+          //---
+          
+          // find data point to be used in Lagrange interpolation (-> k)
+          
+          qeptr = (float **)QE[nPMT];
+          
+          FindLagrange(qeptr,k,wl);
+          
+          // if random > quantum efficiency, reject it
+          
+          qe = Lagrange(qeptr,k,wl) / 100.0;
+          
+          // fprintf(stdout, "%f\n", qe);
+          
+          if ( RandomNumber > qe ) {
+            
+            // read next CPhoton
+            if ( Data_From_STDIN ) 
+              cin.read( flag, SIZE_OF_FLAGS );
+            else
+              inputfile.read ( flag, SIZE_OF_FLAGS );
+            
             // go to beginning of loop
             continue;
@@ -1289 +1196 @@
           
           fnpix[nPMT] += 1.0;
-
+          
 #ifdef __ROOT__ 
           fnslicesum[t_chan]  += 1.0 ; 
           slices[nPMT][t_chan] += 1.0 ; 
 #endif // __ROOT__        
-
+          
 #ifdef __DETAIL_TRIGGER__ 
           //
@@ -1300 +1207 @@
           //
           //
-
+          
           Trigger.Fill( nPMT, ( t - t_ini  ) ) ; 
 #endif // __DETAIL_TRIGGER__ 
           
           // read next CPhoton
-
+          
           if ( Data_From_STDIN ) 
             cin.read( flag, SIZE_OF_FLAGS );
           else
             inputfile.read ( flag, SIZE_OF_FLAGS );
-
+          
         }  // end while, i.e. found end of event
         
         log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
             ncph, ntcph);
-
+      
         // show number of photons
         
         //cout << ncph << " photons read . . . " << endl << flush;
-        
+      
         // skip it ?
         
@@ -1394 +1301 @@
         
         // if we should apply any kind of correction, do it here.
-
+        
         for ( i=0; i<ct_NPixels; ++i ) 
           fnpix[i] *= fCorrection;
-
+        
 #ifdef __DETAIL_TRIGGER__ 
-      //   Trigger.Print() ; 
-      cout << Trigger.Diskriminate() << endl << endl ;
+        //   Trigger.Print() ; 
+        cout << Trigger.Diskriminate() << endl << endl ;
 #endif // __DETAIL_TRIGGER__ 
-
+        
 #ifdef __ROOT__
-
-      //
-      //  Fill the header of this event 
-      //
-      
-      Evt->FillHeader ( (UShort_t) (ntshow + nshow) ,  20 ) ; 
-
-      //  now put out the data of interest
-      //
-      //  1.  -> look for the first slice with signal
-      //
-
-      for ( i=0; i<(2 * SLICES); i++ ) 
-        if ( fnslicesum[i] > 0. )
-          break ; 
-
-      startchan = i ; 
-
-      //
-      //     copy the slices out of the big array 
-      //   
-      //     put the first slice with signal to slice 4
-      //
-      
-      for (i=0; i<ct_NPixels; i++ ) 
-        for ( ii=(startchan-3); ii < (startchan+12); ii++ )  
-          slices2 [i][ii-startchan+3] = slices [i][ii] ; 
-
-
-      // 
-      //  if a pixes has a signal put it to the MRawEvt
-      //
-      
-      for (i=0; i<ct_NPixels; i++ ) {
-        if ( fnpix[i] > 0 ) {
-
-          for ( ii=0; ii < 15; ii++ ) {
-            trans [ii] = slices2[i][ii] ; 
+        
+        //
+        //  Fill the header of this event 
+        //
+        
+        Evt->FillHeader ( (UShort_t) (ntshow + nshow) ,  20 ) ; 
+        
+        //  now put out the data of interest
+        //
+        //  1.  -> look for the first slice with signal
+        //
+        
+        for ( i=0; i<(2 * SLICES); i++ ) 
+          if ( fnslicesum[i] > 0. )
+            break ; 
+        
+        startchan = i ; 
+        
+        //
+        //     copy the slices out of the big array 
+        //   
+        //     put the first slice with signal to slice 4
+        //
+        
+        for (i=0; i<ct_NPixels; i++ ) 
+          for ( ii=(startchan-3); ii < (startchan+12); ii++ )  
+            slices2 [i][ii-startchan+3] = slices [i][ii] ; 
+        
+        
+        // 
+        //  if a pixes has a signal put it to the MRawEvt
+        //
+        
+        for (i=0; i<ct_NPixels; i++ ) {
+          if ( fnpix[i] > 0 ) {
+            
+            for ( ii=0; ii < 15; ii++ ) {
+              trans [ii] = slices2[i][ii] ; 
+            }
+            
+            Evt->FillPixel ( (UShort_t) i , trans ) ; 
+            
           }
-        
-          Evt->FillPixel ( (UShort_t) i , trans ) ; 
-          
         }
-      }
-      
-      //
-      //   
-      //
-      
-      McEvt->Fill( (UShort_t) mcevth.get_primary() , 
-                   mcevth.get_energy(), 
-                   mcevth.get_theta(), 
-                   mcevth.get_phi(), 
-                   mcevth.get_core(),
-                   mcevth.get_coreX(),
-                   mcevth.get_coreY(),
-                   flli,
-                   ulli, ulli, ulli, ulli, ulli ) ; 
-      
-      //
-      //    write it out to the file outfile
-      // 
-
-      EvtTree.Fill() ; 
-
-      //    clear all
-
-      Evt->Clear() ; 
-      McEvt->Clear() ; 
-
+        
+        //
+        //   
+        //
+        
+        McEvt->Fill( (UShort_t) mcevth.get_primary() , 
+                     mcevth.get_energy(), 
+                     mcevth.get_theta(), 
+                     mcevth.get_phi(), 
+                     mcevth.get_core(),
+                     mcevth.get_coreX(),
+                     mcevth.get_coreY(),
+                     flli,
+                     ulli, ulli, ulli, ulli, ulli ) ; 
+        
+        //
+        //    write it out to the file outfile
+        // 
+        
+        EvtTree.Fill() ; 
+        
+        //    clear all
+        
+        Evt->Clear() ; 
+        McEvt->Clear() ; 
+        
 #endif // __ROOT__
         
@@ -1486 +1393 @@
         //--------------------------------------------------
         
-#ifdef __DEBUG__  
-        printf("\n");
-        
-        for ( ici=0; ici<PIX_ARRAY_SIDE; ++ici ) {
-          
-          for ( icj=0; icj<PIX_ARRAY_SIDE; ++icj ) {
-            
-            if ( (int)pixels[ici][icj][PIXNUM] > -1 ) {
-              
-              if ( fnpix[(int)pixels[ici][icj][PIXNUM]] > 0. ) {
-                
-                printf ("@@ %4d %4d %10f %10f %4f (%4d %4d)\n", nshow, 
-                        (int)pixels[ici][icj][PIXNUM], 
-                        pixels[ici][icj][PIXX],
-                        pixels[ici][icj][PIXY],
-                        fnpix[(int)pixels[ici][icj][PIXNUM]], ici, icj);
-                
-              }
-              
-            }
-            
-          }
-          
-        }
-        
-        for (i=0; i<ct_NPixels; ++i) {
-          printf("%d (%d): ", i, npixneig[i]);
-          for (j=0; j<npixneig[i]; ++i) 
-            printf(" %d", pixneig[i][j]);
-          printf("\n");
-        }
-        
-#endif // __DEBUG__
-        
 #ifdef __TRIGGER__
         
@@ -1564 +1437 @@
         //@ If the input parameter "threshold" is 0 we find the maximum 
         //@ trigger threshold this event can pass
-
+        
         for(k=0; ( qThreshold == 0 ? (k <= iMAX_THRESHOLD_PHE) : (k == 1) ); k++){
           
@@ -1666 +1539 @@
                 
                 for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1; ++j ) {
-                  
+                
                   if ( fnpix[pixneig[i][j]] > q0 ) {
                     
@@ -1736 +1609 @@
                     
                     trigger = FALSE;
-                  
+                    
                   }
                   
@@ -1842 +1715 @@
           
           // calculate moments and other things
-
+          
           moments_ptr = moments( anaPixels, fnpixclean, pixary,
-                               plateScale_cm2deg, 0 );
+                                 plateScale_cm2deg, 0 );
           
           charge = moments_ptr->charge ;
@@ -1861 +1734 @@
           asymx  = moments_ptr->asymx  ;
           phiasym= moments_ptr->phi;
-
+          
           lenwid_ptr = lenwid( anaPixels, fnpixclean, pixary,
                                plateScale_cm2deg,
                                ct_PixelWidth_corner_2_corner_half);
-          
-
+          
+          
           // fill the diagnostic Tree
           
@@ -1989 +1862 @@
           image_data[i] =       -1.; i++;
           image_data[i] =       -1.; i++;
-
+          
           // there should be "nvar" variables
-
+          
           if ( i != nvar ) 
             error( SIGNATURE, "Wrong entry length for Ntuple.\n" );
@@ -1999 +1872 @@
           
           // put this information in the data file, 
-
+          
           if ( Write_All_Data ) {
-
+            
             datafile << ntrigger;
             for ( i=0; i<nvar; ++i )
@@ -2016 +1889 @@
           
         } // trigger == FALSE 
-        
+        
 #endif // __TRIGGER__
-
+        
         //!@' @#### Save data.
         //@'
@@ -2027 +1900 @@
         // the output file
         //--------------------------------------------------
-
+        
         //++ 
         // save the image to the file
@@ -2035 +1908 @@
         
         outputfile.write( (char *)&mcevth, mcevth.mysize() ); 
-      
+        
 #ifdef __TRIGGER__
-
+        
         // save the image 
-      
+        
         if ( (trigger == TRUE) || (Write_All_Images == TRUE) )
           outputfile.write( (char *) fnpix, ct_NPixels * sizeof( float ) );
-
+        
 #else
-
+        
         // save the image 
-      
+        
         outputfile.write( (char *) fnpix, ct_NPixels * sizeof( float ) );
-
+        
 #endif // __TRIGGER__
-
+        
         if ( Data_From_STDIN ) 
           cin.read( flag, SIZE_OF_FLAGS );
@@ -2057 +1930 @@
         
       } // end while there is a next event
-
+      
       if( !isA( flag, FLAG_END_OF_RUN   )){
         error( SIGNATURE, "Expected end of run flag, but found: %s\n", flag );
@@ -2101 +1974 @@
             
           }
-        
+          
         } // end if found end of file
       } // end if found end of run
@@ -2110 +1983 @@
     } // end if else found start of run
   } // end big while loop
-
-  //!@' @#### End of program.
-  //@'
-
+  
+//!@' @#### End of program.
+//@'
+  
   //end my version
 
 #ifdef __ROOT__
-      //++
-      // put the Event to the root file
-      //--
-
-      EvtTree.Write() ; 
-      outfile.Write() ;
-      outfile.Close() ; 
+  //++
+  // put the Event to the root file
+  //--
+  
+  EvtTree.Write() ; 
+  outfile.Write() ;
+  outfile.Close() ; 
 
 #endif // __ROOT__
-              
+  
   // close input file
   
@@ -2134 +2007 @@
   log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n", 
        ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
-
+  
   // close files
   
   log( SIGNATURE, "Closing files\n" );
-
+  
   inputfile.close();
   outputfile.close();
@@ -2797 +2670 @@
     /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
 
-    pcam->dpixsizefactor[ipixno] = 1.;
+    pcam->dpixsizefactor[ipixno-1] = 1.;
 
     in = 0;
@@ -3104 +2977 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.3  2000/01/20 18:22:17  petry
+// Found little bug which makes camera crash if it finds a photon
+// of invalid wavelength. This bug is now fixed and the range
+// of valid wavelengths extended to 290 - 800 nm.
+// This is in preparation for the NSB simulation to come.
+// Dirk
+//
 // Revision 1.2  1999/11/19 08:40:42  harald
 // Now it is possible to compile the camera programm under osf1.