Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 340)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 344)
@@ -21,7 +21,7 @@
 //
 // $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,93 +1073,33 @@
 	  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,8 +1111,8 @@
 	  
 	  // 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,88 +1123,20 @@
 	    // 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,4 +1145,39 @@
 	      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,10 +1196,10 @@
 	  
 	  fnpix[nPMT] += 1.0;
-
+	  
 #ifdef __ROOT__ 
 	  fnslicesum[t_chan]  += 1.0 ; 
 	  slices[nPMT][t_chan] += 1.0 ; 
 #endif // __ROOT__ 	  
-
+	  
 #ifdef __DETAIL_TRIGGER__ 
 	  //
@@ -1300,24 +1207,24 @@
 	  //
 	  //
-
+	  
 	  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,84 +1301,84 @@
 	
 	// 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,38 +1393,4 @@
 	//--------------------------------------------------
 	
-#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,5 +1437,5 @@
 	//@ 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,5 +1539,5 @@
 		
 		for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1; ++j ) {
-		  
+		
 		  if ( fnpix[pixneig[i][j]] > q0 ) {
 		    
@@ -1736,5 +1609,5 @@
 		    
 		    trigger = FALSE;
-		  
+		    
 		  }
 		  
@@ -1842,7 +1715,7 @@
 	  
 	  // calculate moments and other things
-
+	  
 	  moments_ptr = moments( anaPixels, fnpixclean, pixary,
-                               plateScale_cm2deg, 0 );
+				 plateScale_cm2deg, 0 );
 	  
 	  charge = moments_ptr->charge ;
@@ -1861,10 +1734,10 @@
 	  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,7 +1862,7 @@
 	  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,7 +1872,7 @@
 	  
 	  // put this information in the data file, 
-
+	  
 	  if ( Write_All_Data ) {
-
+	    
 	    datafile << ntrigger;
 	    for ( i=0; i<nvar; ++i )
@@ -2016,7 +1889,7 @@
 	  
 	} // trigger == FALSE 
-        
+	
 #endif // __TRIGGER__
-
+	
 	//!@' @#### Save data.
 	//@'
@@ -2027,5 +1900,5 @@
 	// the output file
 	//--------------------------------------------------
-
+	
 	//++ 
 	// save the image to the file
@@ -2035,20 +1908,20 @@
 	
 	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,5 +1930,5 @@
 	
       } // 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,5 +1974,5 @@
 	    
 	  }
-	
+	  
 	} // end if found end of file
       } // end if found end of run
@@ -2110,21 +1983,21 @@
     } // 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,9 +2007,9 @@
   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,5 +2670,5 @@
     /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */
 
-    pcam->dpixsizefactor[ipixno] = 1.;
+    pcam->dpixsizefactor[ipixno-1] = 1.;
 
     in = 0;
@@ -3104,4 +2977,11 @@
 //
 // $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.