Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 344)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 365)
@@ -21,7 +21,7 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.4 $
+// $Revision: 1.5 $
 // $Author: petry $ 
-// $Date: 2000-01-25 08:36:23 $
+// $Date: 2000-02-18 17:40:35 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -54,6 +54,4 @@
 
 #include "MDiag.h"
-
-#include "MTrigger.hxx"
 
 #include "MRawEvt.h"
@@ -86,46 +84,4 @@
 #undef  __DEBUG__
 
-// flag for NNT in CT1 camera (default: ON )
-#undef  __CT1_NO_NEIGHBOURS__
-#define __CT1_NO_NEIGHBOURS__
-
-// flag for calculation of NSB (default: ON )
-#undef  __NSB__
-#define __NSB__
-
-// flag for calculation of QE for pixels (default: ON )
-#undef  __QE__
-#define __QE__
-
-
-// flag for implementation of DETAIL_TRIGGER (default: ON )
-//
-//      This is the new implementation of Trigger studies
-//      It relies on a better simulation of the time stucture 
-//      of the PhotoMultiplier. For more details see the 
-//      documentation of the --> class MTrigger <-- 
-#define __DETAIL_TRIGGER__
-#undef  __DETAIL_TRIGGER__
-
-// flag for implementation of TRIGGER (default: ON )
-#undef  __TRIGGER__
-#define __TRIGGER__
-
-// flag for implementation of Tail-Cut (default: ON )
-#undef  __TAILCUT__
-#define __TAILCUT__
-
-// flag for calculation of islands stat. (default: ON )
-#define __ISLANDS__
-#undef  __ISLANDS__
-
-// flag for calculation of image parameters (default: ON )
-#undef  __MOMENTS__
-#define __MOMENTS__
-
-// flag for ROOT  (default: ON )
-#undef  __ROOT__
-#define __ROOT__
-
 //!@}
 
@@ -259,6 +215,4 @@
 // and data stored on them with information about the pixels
 
-//@: table for IJ sys.
-static float pixels[PIX_ARRAY_SIDE][PIX_ARRAY_SIDE][4];   
 
 //@: coordinates x,y for each pixel
@@ -276,7 +230,4 @@
 //@: contents of the pixels (ph.e.) after cleanning
 static float *fnpixclean; 
-
-//@: contents of the sum of all ph.e. in one timeslice of 1 ns
-static float *fnslicesum ;
 
  
@@ -364,30 +315,4 @@
   @"*/
 
-//!@{
-// table of random numbers
-
-// (unused)
-//static double RandomNumbers[500];  
-//!@}
-
-/*!@"
-
-  The following is a variable to count the number of Cphotons
-  in the different steps of the simulation. 
-  The definition is as follows:
-  @[
-  \mbox{CountCphotons}[ \mbox{FILTER} ] \equiv
-  \mbox{\it Number of photons after the filter} \mbox{FILTER}
-  @]
-  The filters are defined and can be found in the file |camera.h|.
-
-  @"*/
-
-//!@{
-// vector to count photons at any given step of the simulation
-
-static int CountCphotons[10];  
-//!@}
-
 /*!@"
 
@@ -430,5 +355,5 @@
 
   char inname[256];           //@< input file name
-  char outname[256];          //@< output file name
+  char starfieldname[256];    //@< starfield input file name
   char datname[256];          //@< data (ASCII) output file name
   char diagname[256];         //@< diagnistic output file (ROOT format)
@@ -437,14 +362,16 @@
   char parname[256];          //@< parameters file name
 
-  char sign[20];              //@< initialize sign
-
   char flag[SIZE_OF_FLAGS + 1];  //@< flags in the .rfl file
 
-  ifstream inputfile;         //@< stream for the input file
-  ofstream outputfile;        //@< stream for the output file
+  FILE *inputfile;            //@< stream for the input file
   ofstream datafile;          //@< stream for the data file
 
   MCEventHeader mcevth;       //@< Event Header class (MC)
-  MCCphoton cphoton;          //@< Cherenkov Photon class (MC)
+
+  Photoelectron *photoe = NULL; //@< array of the photoelectrons of one event
+  int inumphe;                //@< number of photoelectrons in an event
+
+  float arrtmin_ns;           //@ arrival time of the first photoelectron
+  float arrtmax_ns;           //@ arrival time of the last photoelectron
 
   float thetaCT, phiCT;       //@< parameters of a given shower
@@ -458,26 +385,25 @@
   int nshow=0;                //@< partial number of shower in a given run
   int ntshow=0;               //@< total number of showers
-  int ncph=0;                 //@< partial number of shower in a given run
-  int ntcph=0;                //@< total number of showers
-
-  int i, ii, j, k;            //@< simple counters
-
-  float t_ini;                //@< time of the first Cphoton read in
-  float t;                    //@< time for a single photon
-  int   t_chan ;              //@< the bin (channel) in time of a single photon
-  
-  int   startchan ;           //@< the first bin with entries in the time slices
-
-  float cx, cy, ci, cj;       //@< coordinates in the XY and IJ systems
-  int ici, icj, iici, iicj;   //@< coordinates in the IJ (integers)
-
-  int nPMT;                   //@< number of pixel
-
-  float wl, last_wl;          //@< wavelength of the photon
-  float qe;                   //@< quantum efficiency
-  float **qeptr;
+  int ncph=0;                 //@< partial number of photons in a given run
+  int ntcph=0;                //@< total number of photons
+
+  int i, j, k;                //@< simple counters
 
   int simulateNSB;            //@< Will we simulate NSB?
-  float meanNSB;              //@< NSB mean value (per pixel)
+  float meanNSB;              //@< diffuse NSB mean value (phe per ns per central pixel)
+  float diffnsb_phepns[iMAXNUMPIX];  //@< diffuse NSB values for each pixel derived 
+                                     //@< from meanNSB
+  float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS];    //@< non-diffuse nsb 
+                                                     //@< photoelectron rates 
+  float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
+    EXTWAVEBAND1,
+    EXTWAVEBAND2,
+    EXTWAVEBAND3,
+    EXTWAVEBAND4,
+    EXTWAVEBAND5
+  };                          
+  float baseline_mv[iMAXNUMPIX]; //@< The baseline (mV) caused by the NSB; to be subtracted
+                                 //@< in order to simulate the preamps' AC coupling 
+
   float qThreshold;           //@< Threshold value
   float qTailCut;             //@< Tail Cut value
@@ -488,14 +414,7 @@
   float fCorrection;          //@< Factor to apply to pixel values (def. 1.)
 
-  float q0;                   //@< trigger threshold ( intermediate variable )
-  float maxcharge;            //@< maximum charge in pixels
-  int noverq0, novq0;         //@< number of pixels above threshold
-  int ngrpq0, mxgrp;          //@< number of pixels in a group
-
   int trigger;                //@< trigger flag 
   int itrigger;               //@< index of pixel fired
   int ntrigger = 0;           //@< number of triggers in the whole file
-  unsigned char triggerBits;  //@< byte for trigger condition check (MAGIC)
-  int bit;                    //@< intermediate variable
 
   float plateScale_cm2deg;    //@< plate scale (deg/cm)
@@ -506,6 +425,4 @@
   int still_in_loop = FALSE;
 
-  char    Signature[20];
-
   float *image_data;
   int nvar, hidt;
@@ -513,9 +430,4 @@
   struct camera cam; // structure holding the camera definition
 
-#ifdef __DETAIL_TRIGGER__
-
-  MTrigger  Trigger ;         //@< A instance of the Class MTrigger 
-
-#endif // __DETAIL_TRIGGER__ 
 
   //!@' @#### Definition of variables for |getopt()|.
@@ -592,5 +504,5 @@
 
   strcpy( inname, get_input_filename() );
-  strcpy( outname, get_output_filename() );
+  strcpy( starfieldname, get_starfield_filename() );
   strcpy( datname, get_data_filename() );
   strcpy( diagname, get_diag_filename() );
@@ -615,9 +527,10 @@
       "Filenames",
       "In", inname, 
-      "Out", outname, 
+      "Stars", starfieldname,
+      "CT", ct_filename, 
       "Data", datname,
       "Diag", diagname,
-      "ROOT",  rootname, 
-      "CT", ct_filename);
+      "ROOT",  rootname 
+      );
 
   
@@ -669,8 +582,12 @@
   read_ct_file();
 
-  // read pixels data
+  // read camera setup
 
   read_pixels(&cam);
 
+  // allocate memory for the photoelectrons
+
+  photoe = new Photoelectron[iMAXNUMPHE];
+
   // initialise ROOT
 
@@ -682,5 +599,4 @@
 
   hfile = new TFile( diagname,"RECREATE", "MAGIC Telescope MC diagnostic data");
-  
   
   // Create the ROOT Tree for the diagnostic data
@@ -698,5 +614,5 @@
   tree->Branch("event", "MDiagEventobject", &event, bsize, split);
 
-#ifdef __ROOT__
+  // Prepare the raw data output
 
   MRawEvt *Evt   = new MRawEvt() ; 
@@ -706,9 +622,7 @@
   //
 
-  TFile outfile ( rootname , "RECREATE" ) ; 
-
-  //
+  TFile outfile ( rootname , "RECREATE" ); 
+
   //      create a Tree for the Event data stream 
-  //
 
   TTree EvtTree("EvtTree","Events of Run");
@@ -722,9 +636,4 @@
   		 &McEvt, bsize, split);
 
-
-  float  flli = 0. ; 
-  unsigned short ulli = 0 ; 
-
-#endif // __ROOT__
 
   // for safety and for dimensioning image_data: count the elements in the 
@@ -787,59 +696,72 @@
   anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;
 
-  // open input file if we DO read data from a file
-
-  if (! Data_From_STDIN) {  
-    log( SIGNATURE, "Openning input \"rfl\" file %s\n", inname );
-    inputfile.open( inname );
-    if ( inputfile.bad() ) 
+  // prepare the NSB simulation
+
+  if( simulateNSB ){
+
+    //
+    // Calculate the non-diffuse NSB photoelectron rates
+    //
+    
+    k = produce_nsbrates( starfieldname,
+			  &cam,
+			  photoe, // only a dummy here
+			  nsbrate_phepns );
+    if (k != 0){
+      cout << "Error when reading starfield... \nExiting.\n";
+      exit(1);
+    }
+  
+//     for(i=0; i<cam.inumpixels; i++){
+//       cout << i;
+//       for(j=0; j<iNUMWAVEBANDS; j++){
+// 	   cout << " " << j << " " << nsbrate_phepns[i][j];
+//       }
+//       cout << "\n";
+//     }
+
+    // calculate diffuse rate correcting for the pixel size
+
+    for(i=0; i<cam.inumpixels; i++){
+      diffnsb_phepns[i] = meanNSB * 
+	cam.dpixsizefactor[i] * cam.dpixsizefactor[i];
+    }
+
+  }
+
+  //
+  // Read the reflector file with the Cherenkov data
+  //			
+
+  // select input file 
+
+  if ( Data_From_STDIN ) {
+
+    inputfile = stdin;
+
+  }
+  else{
+
+    log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname );
+    inputfile = fopen( inname, "r" );
+    if ( inputfile == NULL ) 
       error( SIGNATURE, "Cannot open input file: %s\n", inname );
+
   }
   
   // get signature, and check it
-  // NOTE: this part repeats further down in the code;
-  // if you change something here you probably want to change it 
-  // there as well
-
-  strcpy(Signature, REFL_SIGNATURE);
-
-  strcpy(sign, Signature);
-
-  if ( Data_From_STDIN ) 
-    cin.read( (char *)sign, strlen(Signature));
-  else
-    inputfile.read( (char *)sign, strlen(Signature));
-
-  if (strcmp(sign, Signature) != 0) {
-    cerr << "ERROR: Signature of .rfl file is not correct\n";
-    cerr << '"' << sign << '"' << '\n';
-    cerr << "should be: " << Signature << '\n';
+
+  if(check_reflector_file( inputfile )==FALSE){
     exit(1);
   }
 
-  if ( Data_From_STDIN ) 
-    cin.read( (char *)sign, 1);
-  else
-    inputfile.read( (char *)sign, 1);
-
-  // open output file
-  
-  log( SIGNATURE, "Openning output \"phe\" file %s\n", outname );
-  outputfile.open( outname );
-  
-  if ( outputfile.bad() ) 
-    error( SIGNATURE, "Cannot open output file: %s\n", outname );
-  
   // open data file
   
-  log( SIGNATURE, "Openning data \"dat\" file %s\n", datname );
+  log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
   datafile.open( datname );
   
-  if ( outputfile.bad() ) 
-    error( SIGNATURE, "Cannot open output file: %s\n", outname );
-  
-  // write signature
-
-  outputfile.write( SIGNATURE, sizeof(SIGNATURE) );
-
+  if ( datafile.bad() ) 
+    error( SIGNATURE, "Cannot open data file: %s\n", datname );
+  
   // initializes flag
   
@@ -851,28 +773,18 @@
   fnpixclean = new float [ ct_NPixels ];
 
-#ifdef __ROOT__
-  
-  fnslicesum = new float [ (2 * SLICES) ] ; 
-  
-  float slices   [CAMERA_PIXELS][ (2 * SLICES) ] ; 
-  float slices2  [CAMERA_PIXELS][ SLICES ] ; 
-
-  float trans    [ SLICES ] ; 
-#endif // __ROOT__ 
-
-  
   moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 );
+
+  // initialize baseline 
+
+  for(i=0; i<cam.inumpixels; i++){
+    baseline_mv[i] = 0.;
+  }
         
   //!@' @#### Main loop.
   //@'
 
-  //begin my version
-                                           
   // get flag
     
-  if ( Data_From_STDIN ) 
-    cin.read( flag, SIZE_OF_FLAGS );
-  else
-    inputfile.read ( flag, SIZE_OF_FLAGS );
+  fread( flag, SIZE_OF_FLAGS, 1, inputfile );
 
   // loop over the file
@@ -881,5 +793,5 @@
 
   while (
-         ((! Data_From_STDIN) && (! inputfile.eof()))
+         ((! Data_From_STDIN) && ( !feof(inputfile) ))
          ||
          (Data_From_STDIN && still_in_loop)
@@ -891,41 +803,10 @@
     }
     else { // found start of run
+
       nshow=0;
-      // read next flag
-
-      if ( Data_From_STDIN ) 
-	cin.read( flag, SIZE_OF_FLAGS );
-      else
-	inputfile.read ( flag, SIZE_OF_FLAGS );
+
+      fread( flag, SIZE_OF_FLAGS, 1, inputfile );
 
       while( isA( flag, FLAG_START_OF_EVENT   )){ // while there is a next event
-	/*!@'
-	  
-	  For the case  |FLAG\_START\_OF\_EVENT|,
-	  we read each Cherenkov photon, and follow these steps:
-	  
-	  @enumerate
-	  
-	  @- Transform XY-coordinates to IJ-coordinates.
-	  
-	  @- With this, we obtain the pixel where the photon hits.
-	  
-	  @- Use the wavelength $\lambda$ and the table of QE, and
-	  calculate the estimated (third order interpolated) quantum
-	  efficiency for that photon. The photon can be rejected.
-	  
-	  @- If accepted, then add to the pixel.
-	  
-	  @endenumerate
-	  
-	  In principle, we should stop here, and use another program to
-	  'smear' the image, to add the Night Sky Background, and to
-	  simulate the trigger logic, but we will make this program
-	  quick and dirty, and include all here.
-	  
-	  If we are reading PHE files, we jump to the point where the
-	  pixelization process already has finished.
-	  
-	*/
 	
 	++nshow;
@@ -934,8 +815,5 @@
 	// get MCEventHeader
 	
-	if ( Data_From_STDIN ) 
-	  cin.read( (char*)&mcevth, mcevth.mysize() );
-	else
-	  mcevth.read( inputfile );
+	fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
 	
 	// calculate core distance and impact parameter
@@ -1023,208 +901,28 @@
 	}
 
-	// clear camera
-	
-	for ( i=0; i<ct_NPixels; ++i ){
- 
-	  fnpix[i] = 0.0;
-#ifdef __ROOT__	
-	  for ( ii=0; ii<(2 * SLICES); ii++ ) { 
-	    slices [i][ii] = 0 ; 
-	  }
-#endif // __ROOT__ 
+	// read the photons and produce the photoelectrons
+
+	k = produce_phes( inputfile,
+			  &cam,
+			  WAVEBANDBOUND1,
+			  WAVEBANDBOUND6,
+			  photoe,   // will be changed by the function!
+			  &inumphe, // important for later: the size of photoe[]
+			  fnpix,    // will be changed by the function!
+			  &ncph,    // will be changed by the function!
+			  &arrtmin_ns, // will be changed by the function!
+			  &arrtmax_ns // will be changed by the function!
+			  );
+
+	if( k != 0 ){ // non-zero returnvalue means error
+	  cout << "Exiting.\n";
+	  exit(1);
 	}
-
-	ntcph +=ncph;
-	ncph = 0;
-
-#ifdef __DETAIL_TRIGGER__ 
-	//
-	//   clear Trigger 
-	//
-      
-	Trigger.Reset() ; 
-#endif // __DETAIL_TRIGGER__ 
-       
-	//- - - - - - - - - - - - - - - - - - - - - - - - - 
-	// read photons and "map" them into the pixels
-	//--------------------------------------------------      
-	
-	// initialize CPhoton
-	
-	cphoton.fill(0., 0., 0., 0., 0., 0., 0., 0.);
-	
-	// read the photons data
-	
-	if ( Data_From_STDIN ) 
-	  cin.read( flag, SIZE_OF_FLAGS );
-	else
-	  inputfile.read ( flag, SIZE_OF_FLAGS );
-	 
-	// loop over the photons
-
-	t_ini = -99999;
-	
-	while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
-	  
-	  memcpy( (char*)&cphoton, flag, SIZE_OF_FLAGS );
-
-	  if ( Data_From_STDIN ) 
-	    cin.read( ((char*)&cphoton)+SIZE_OF_FLAGS, cphoton.mysize()-SIZE_OF_FLAGS );
-	  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);
-	  }
-	  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;
-	  } 
-	    
-	  // Pixelization
-	    
-	  cx = cphoton.get_x();
-	  cy = cphoton.get_y(); 
-	  
-	  // get wavelength
-	  
-	  last_wl = wl;
-	  wl = cphoton.get_wl();
-	  
-	  // 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 ) 
-	      cin.read( flag, SIZE_OF_FLAGS );
-	    else
-	      inputfile.read ( flag, SIZE_OF_FLAGS );
-	    
-	    // go to beginning of loop, the photon is lost
-	    continue;
-	    
-	  }
-
-	  // cout << "@#1 " << nshow << ' ' << cx << ' ' << cy << endl;
-	  
-	  nPMT = -1;
-
-	  for(i=0; 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 ) {
-	    
-	    // 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;
-	    
-	  }
-	  
-#endif // __QE__
-	  
-	  //+++
-	  // Cphoton is accepted
-	  //---
-	  
-	  // increase the number of Cphs. in the PMT, i.e.,
-	  // increase in one unit the counter of the photons
-	  // stored in the pixel nPMT
-	  
-	  fnpix[nPMT] += 1.0;
-	  
-#ifdef __ROOT__ 
-	  fnslicesum[t_chan]  += 1.0 ; 
-	  slices[nPMT][t_chan] += 1.0 ; 
-#endif // __ROOT__ 	  
-	  
-#ifdef __DETAIL_TRIGGER__ 
-	  //
-	  //  fill the Trigger class with this phe
-	  //
-	  //
-	  
-	  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;
-      
+
+	ntcph += ncph;
+
 	// skip it ?
 	
@@ -1245,32 +943,5 @@
 	}
 	
-	/*!@'
-	  
-	  After reading all the Cherenkov photons for a given event,
-	  we have in the table of number of photons for each pixel
-	  only the 'raw' amount of Cherenkov photons @$n_p@$. Now, we
-	  should take this number as the mean value of the
-	  distribution of photons in that pixel @$p@$, following a
-	  Poisson distribution.
-	  
-	  @[ n_p \equiv \mu_p @]
-	  
-	  and with this number the amount of light coming from the
-	  shower is calculated @$\hat{n}_p@$.
-	  
-	  Then, we calculate the amount of Night Sky Background we
-	  must introduce in that pixel @$p@$. We calculate this using
-	  again a Poisson distribution with mean @$\mu_\mathrm{NSB}@$
-	  (defined in the |camera.h| file). The value of
-	  @$\mu_\mathrm{NSB}@$ is obtained from measurements. With this
-	  value, the amount of photons @$\hat{n}_\mathrm{NSB}@$ coming
-	  from the Night Sky Background is calculated.
-	  
-	  Finally, the amount of photons for that pixels is:
-	  @[ \hat{n}_p^\mathrm{final} = \hat{n}_p + \hat{n}_\mathrm{NSB} @]
-	  
-	*/
-	
-	// after reading all the photons, our camera is filled
+	// energy cut
 	
 	if ( Select_Energy ) {
@@ -1281,375 +952,47 @@
 	  }
 	}
-	
-#ifdef __NSB__
-	
-	//!@' @#### NSB (Night Sky Background) simulation.
-	//@'
-	
-	//+++
+
 	// NSB simulation
-	//---
-	
-	// add NSB "noise"
-	// TO DO: make meanNSB an array and read the contents from a file!
-	
-	if ( simulateNSB )
-	  for ( i=0; i<ct_NPixels; ++i ) 
-	    fnpix[i] += (float)ignpoi( meanNSB );
-	
-#endif // __NSB__
-	
-	// 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 ;
-#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] ; 
-	    }
-	    
-	    Evt->FillPixel ( (UShort_t) i , trans ) ; 
-	    
+
+	if(simulateNSB){
+
+	  k = produce_nsb_phes( &arrtmin_ns, // will be changed by the function!
+				&arrtmax_ns, // will be changed by the function!
+				thetaCT,
+				&cam,
+				nsbrate_phepns,
+				diffnsb_phepns,
+				ext,
+				fnpix,   // will be changed by the function!
+				photoe,  // will be changed by the function!
+				&inumphe, // important for later: the size of photoe[]
+				baseline_mv // will be generated by the function
+				);
+
+	  if( k != 0 ){ // non-zero returnvalue means error
+	    cout << "Exiting.\n";
+	    exit(1);
 	  }
-	}
-	
-	//
-	//   
-	//
-	
-	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__
-	
-	//++++++++++++++++++++++++++++++++++++++++++++++++++
-	// at this point we have a camera full of
-	// ph.e.s
-	// we should first apply the trigger condition,
-	// and if there's trigger, then clean the image,
-	// calculate the islands statistics and the
-	// other parameters of the image (Hillas' parameters
-	// and so on).
-	//--------------------------------------------------
-	
-#ifdef __TRIGGER__
-	
-	/*!@'
-	  
-	  @#### Trigger logic simulation.
-	  
-	  In the following block we look at the pixel contents, looking
-	  for pixels fulfilling the trigger condition. This condition,
-	  in this current version of the program, is the following:
-	  
-	  @itemize
-	  
-	  @- |CT1|: Two neighbour pixels with charge above the threshold
-	  @$q_0@$. For the old CT1 data, however, the trigger condition
-	  was 'any two pixels with charge above the threshold @$q_0@$'.
-	  
-	  @- |MAGIC|: A 'closed-packet' of four neighbour pixels, each
-	  of them with charge above the threshold @$q_0@$.
-	  
-	  @enditemize
-	  
-	  In the following figure you can find a sort of description 
-	  about the meanning of 'closed-packet'.
-	  
-	  @F
-	  
-	  \begin{figure}[htbp]
-	  \begin{center}
-	  \includegraphics{closepck.eps}
-	  \caption{Meanning of the expression ``{\it close-packet}''}
-	  \label{fig:closepacket}
-	  \end{center}
-	  \end{figure}
-	  
-	  @F
-	  
-	*/
-	
-	//++ 
-	// TRIGGER Condition
-	//--
-	
-	//@ 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++){
-	  
-	  // is there trigger?
-	  
-	  noverq0 = 0;
-	  q0 = ( qThreshold == 0. ? (float) k : qThreshold );
-	  trigger = FALSE;
-	  mxgrp = 0;
-	  maxcharge = 0.0;
-	  
-	  // Warning! NOT all the camera is able to give trigger
-	  // only up to 'degTrigger' degrees 
-	  
-	  for ( i=0 ; (i<ct_NCentralPixels) && (trigger==FALSE) ; ++i ) {
-	    
-	    // calculate absolute maximum
-	    
-	    maxcharge = MAX(fnpix[i],maxcharge); 
-	    
-	    // is this pixel above threshold ?
-	    
-	    if ( fnpix[i] <= q0 )   
-	      continue;
-	    
-	    // it is: increment the number of pixels above threshold
-	    
-	    ++noverq0;
-	    
-	    // if the trigger already fired, just count the pixels 
-	    // above threshold 
-	    
-	    if ( trigger == TRUE )    
-	      continue;
-	    
-	    // is this pixel inside the trigger zone in the camera ?
-	    
-	    if ( (sqrt(SQR(pixary[i][0]) + 
-		       SQR(pixary[i][1]))*plateScale_cm2deg) > degTriggerZone) 
-	      continue;
-	    
-	    // 'ngrpq0' is the number of neighbours of pixel i with q > q0
-	    
-	    ngrpq0 = 0;
-	    
-	    // look at each pixel in the neighborhood, and count
-	    // those above threshold q0
-	    
-	    for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1 ; ++j )
-	      if ( fnpix[pixneig[i][j]] > q0 )
-		++ngrpq0;
-	    
-	    // check whether we have trigger
-	    
-	    if ( ct_Type == 0 ) {
-	      
-	      //++ >>>>> CT1 <<<<<
-	      
-#ifdef __CT1_NO_NEIGHBOURS__
-	      
-	      if ( noverq0 > 1 )
-		trigger = TRUE;
-	      
-#else
-	      
-	      if ( ngrpq0 > 0 )
-		trigger = TRUE;
-	      
-#endif
-	      
-	      //-- >>>>> CT1 <<<<<
-	      
-	    } else {
-	      
-	      //++ >>>>> MAGIC <<<<<
-	      
-	      // (at least 4 packed with at least q0 phes)
-	      
-	      // there are 3 cases
-	      // 1. zero, one or two neighbours have enough charge: no trigger
-	      // 2. five or six neighbours with enough charge: trigger? sure!!
-	      // 3. three or four neighbours with q > q0 : we must look
-	      //    for 'closeness'.
-	      
-	      switch ( ngrpq0 ) {
-		
-	      case 0:
-	      case 1:
-	      case 2:
-		
-		trigger = FALSE;
-		break;
-		
-	      case 3:
-	      case 4:
-		
-		// if reaches this line, it means three or four neighbours
-		// around the central pixel
-		
-		triggerBits = 1;
-		
-		for ( j=0 ; j<npixneig[i] && pixneig[i][j]>-1; ++j ) {
-		
-		  if ( fnpix[pixneig[i][j]] > q0 ) {
-		    
-		    if ( pixary[pixneig[i][j]][0] > pixary[i][0] ) {
-		      
-		      if ( nint(pixary[pixneig[i][j]][1]*10.0) >
-			   nint(pixary[i][1]*10.0) )
-			bit = 2;
-		      else if ( nint(pixary[pixneig[i][j]][1]*10.0) <
-				nint(pixary[i][1]*10.0) )
-			bit = 6;
-		      else 
-			bit = 1;
-		      
-		    } else {
-		      
-		      if ( nint(pixary[pixneig[i][j]][1]*10.0) >
-			   nint(pixary[i][1]*10.0) )
-			bit = 3;
-		      else if ( nint(pixary[pixneig[i][j]][1]*10.0) <
-				nint(pixary[i][1]*10.0) )
-			bit = 5;
-		      else 
-			bit = 4;
-		      
-		    }
-		    
-		    triggerBits |= (1<<bit);
-		    
-		  }
-		  
-		}
-		
-		if ( ngrpq0 == 3 ) {  // 4-fold trigger
-		  
-		  switch ( triggerBits ) {
-		    
-		  case 0x0f:          // 0 000111 1
-		  case 0x1d:          // 0 001110 1
-		  case 0x39:          // 0 011100 1
-		  case 0x71:          // 0 111000 1
-		  case 0x63:          // 0 110001 1
-		  case 0x47:          // 0 100011 1
-		    
-		    trigger = TRUE;
-		    break;
-		    
-		  default:
-		    
-		    trigger = FALSE;
-		    
-		  }
-		  
-		} else {              // 4-fold trigger
-		  
-		  switch ( triggerBits ) {
-		    
-		  case 0x1f:          // 0 001111 1
-		  case 0x3d:          // 0 011110 1
-		  case 0x79:          // 0 111100 1
-		  case 0x73:          // 0 111001 1
-		  case 0x67:          // 0 110011 1
-		  case 0x4f:          // 0 100111 1
-		    
-		    trigger = TRUE;
-		    break;
-		    
-		  default:
-		    
-		    trigger = FALSE;
-		    
-		  }
-		  
-		}
-		
-		mxgrp = MAX(ngrpq0,mxgrp);
-		
-		break;
-		
-	      case 5:
-	      case 6:
-		
-		trigger = TRUE;
-		break;
-		
-	      default:
-		
-		trigger = FALSE;
-		error( SIGNATURE, "Number of neighbours > 6 !!! Exiting.\n\n");
-		break;
-		
-	      } // switch (ngrpq0)
-	      
-	    } // ct_Type
-	    
-	  } // for each pixel i
-	  
-	  if ( trigger == FALSE ) {
-	    break;
-	  } // end if
-	  
-	} //end  for each threshold
-	maxtrigthr_phe = (float) k-1;  // i.e. maxtrigthr_phe < 0. means that
-	// the event doesn't even pass threshold 0.
-	//  maxtrigthr_phe >= 0 means, the event passes some threshold
-	//  or (in case the input parameter "threshold" was > 0), the event
-	//  passes the threshold given by the input parameter. 
-	if ( maxtrigthr_phe >= 0. ) {
-	  trigger = TRUE;
-	}
-	
-	
-	novq0 = noverq0;
+
+	}// end if(simulateNSB) ...
+
+
+//	cout << arrtmin_ns << " " << arrtmax_ns << "\n";
+//	for(i=0; i<cam.inumpixels; i++){
+//        cout << i << " " << baseline_mv[i] <<"\n";
+//      }
+  
+
+	cout << "Total number of phes: " << inumphe << "\n";
+
+	// TRIGGER HERE
+
+	trigger = FALSE;
+
+	if( TRUE ) trigger = TRUE; // put your trigger function here
+
+// 	for( i=0; i<inumphe; i++){
+// 	  cout << "phe " << photoe[i].ipixnum << " " << photoe[i].iarrtime_ns << "\n";
+// 	}
 	
 	if ( trigger == TRUE ) {
@@ -1660,33 +1003,4 @@
 	  memcpy( fnpixclean, fnpix, sizeof(float) * ct_NPixels );
 	  
-#ifdef __TAILCUT__
-	  
-	  //!@' @#### Tail-cut condition.
-	  //@'
-	  
-	  //++
-	  // tail-cut 
-	  //--
-	  
-	  // Tail-Cut = 0   : No Tail-Cut
-	  // Tail-Cut > 0   : Make Tail-Cut
-	  // Tail-Cut < 0   : Make Tail-Cut with t_0 = Sqrt[ maximum ]
-	  
-	  if (qTailCut > 0.0) {
-	    
-	    for ( i=0; i<ct_NPixels; ++i ) 
-	      if ( fnpix[i] < qTailCut ) 
-		fnpixclean[i] = 0.0;
-	    
-	  } else if (qTailCut < 0.0) {
-	    
-	    maxcharge = sqrt(maxcharge);
-	    for ( i=0; i<ct_NPixels; ++i ) 
-	      if ( fnpix[i] < maxcharge ) 
-		fnpixclean[i] = 0.0;
-	    
-	  }
-	  
-#endif // __TAILCUT__
 	  
 #ifdef __ISLANDS__
@@ -1705,5 +1019,4 @@
 #endif // __ISLANDS__
 	  
-#ifdef __MOMENTS__
 	  
 	  //!@' @#### Calculation of parameters of the image.
@@ -1797,8 +1110,4 @@
 	  }
 	  
-	  
-#endif // __MOMENTS__
-	  
-	  
 	  // revert the fnpixclean matrix into fnpix
 	  // (now we do this, but maybe in a future we want to
@@ -1889,7 +1198,5 @@
 	  
 	} // trigger == FALSE 
-	
-#endif // __TRIGGER__
-	
+		
 	//!@' @#### Save data.
 	//@'
@@ -1904,29 +1211,11 @@
 	// save the image to the file
 	//--
-	
-	// write MCEventHeader to output file
-	
-	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 );
-	else
-	  inputfile.read ( flag, SIZE_OF_FLAGS );
-	
+
+
+
+	// look for the next event
+	
+	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
+
       } // end while there is a next event
       
@@ -1938,8 +1227,5 @@
 	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
 	
-	if ( Data_From_STDIN ) 
-	  cin.read( flag, SIZE_OF_FLAGS );
-	else
-	  inputfile.read ( flag, SIZE_OF_FLAGS );
+	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
 	
 	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
@@ -1947,47 +1233,23 @@
 	  still_in_loop  = FALSE;
 	  
-	  if ((! Data_From_STDIN) && (! inputfile.eof())){
+	  if ((! Data_From_STDIN) && ( !feof(inputfile) )){
 	    
 	    // we have concatenated input files.
 	    // get signature of the next part and check it.
-	    // NOTE: this part repeats further up in the code;
-	    // if you change something here you probably want to change it 
-	    // there as well
-	    
-	    strcpy(Signature, REFL_SIGNATURE);
-	    
-	    strcpy(sign, Signature);
-	    
-	    inputfile.read( (char *)sign, strlen(Signature));
-	    
-	    if (strcmp(sign, Signature) != 0) {
-	      cerr << "ERROR: Signature of .rfl file is not correct\n";
-	      cerr << '"' << sign << '"' << '\n';
-	      cerr << "should be: " << Signature << '\n';
+
+	    if(check_reflector_file( inputfile )==FALSE){
 	      exit(1);
 	    }
 	    
-	    if ( Data_From_STDIN ) 
-	      cin.read( (char *)sign, 1);
-	    else
-	      inputfile.read( (char *)sign, 1);
-	    
 	  }
 	  
 	} // end if found end of file
       } // end if found end of run
-      if ( Data_From_STDIN ) 
-	cin.read( flag, SIZE_OF_FLAGS );
-      else
-	inputfile.read ( flag, SIZE_OF_FLAGS );
+
+      fread( flag, SIZE_OF_FLAGS, 1, inputfile );
+
     } // end if else found start of run
   } // end big while loop
   
-//!@' @#### End of program.
-//@'
-  
-  //end my version
-
-#ifdef __ROOT__
   //++
   // put the Event to the root file
@@ -1998,9 +1260,7 @@
   outfile.Close() ; 
 
-#endif // __ROOT__
   
   // close input file
   
-  ntcph += ncph;
   log( SIGNATURE, "%d event(s), with a total of %d C.photons\n", 
        ntshow, ntcph );
@@ -2012,6 +1272,7 @@
   log( SIGNATURE, "Closing files\n" );
   
-  inputfile.close();
-  outputfile.close();
+  if( ! Data_From_STDIN ){
+    fclose( inputfile );
+  }
   datafile.close();
 
@@ -2019,11 +1280,4 @@
 
   hfile->Close();
-
-#ifdef __DETAIL_TRIGGER__ 
-  // Output of Trigger statistics
-  //
-
-  Trigger.PrintStat() ; 
-#endif // __DETAIL_TRIGGER__ 
 
   // program finished
@@ -2138,9 +1392,9 @@
 
   //  Display the name of the function that called error
-  fprintf(stderr, "ERROR in %s: ", funct);
+  fprintf(stdout, "ERROR in %s: ", funct);
 
   // Display the remainder of the message
   va_start(args, fmt);
-  vfprintf(stderr, fmt, args);
+  vfprintf(stdout, fmt, args);
   va_end(args);
 
@@ -2425,5 +1679,5 @@
   ifstream qefile;
   char line[LINE_MAX_LENGTH];
-  int n, i, j, k;
+  int n, i, j, icount;
   float qe;
 
@@ -2438,8 +1692,4 @@
 
   // initialize pixel numbers
-
-  for ( i=0; i<PIX_ARRAY_SIDE; ++i ) 
-    for ( j=0; j<PIX_ARRAY_SIDE; ++j ) 
-      pixels[i][j][PIXNUM] = -1;
 
   pixary = new float* [2*ct_NCentralPixels];
@@ -2462,19 +1712,4 @@
   igen_pixel_coordinates(pcam);
 
-  // transfer coordinates to the working arrays for
-  // the central pixels
-
-  for(k=0; k<ct_NCentralPixels; k++){
-
-    i = (int) pcam->di[k];
-    j = (int) pcam->dj[k];
-
-    pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXNUM] = k;
-    pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXX] = pcam->dxc[k];
-    pixels[i+PIX_ARRAY_HALF_SIDE][j+PIX_ARRAY_HALF_SIDE][PIXY] = pcam->dyc[k];
-   
-    pixary[k][0] = pcam->dxc[k];
-    pixary[k][1] = pcam->dyc[k];
-  }
 
   // calculate tables of neighbours
@@ -2519,5 +1754,5 @@
   // try to open the file
 
-  log("read_pixels", "Openning the file \"%s\" . . .\n", QE_FILE);
+  log("read_pixels", "Opening the file \"%s\" . . .\n", QE_FILE);
   
   qefile.open( QE_FILE );
@@ -2533,4 +1768,5 @@
 
   i=-1;
+  icount = 0;
 
   while ( ! qefile.eof() ) {          
@@ -2577,5 +1813,6 @@
     // get the values (num-pixel, num-datapoint, QE-value)
     
-    sscanf(line, "%d %d %f", &i, &j, &qe);
+    if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 ) 
+      break;
 
     if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
@@ -2585,9 +1822,30 @@
     }
 
+    if ( i > ct_NPixels) break;
+
+    icount++;
+
   }
 
+  if(icount/pointsQE < ct_NPixels){
+    error( "read_pixels", "The quantum efficiency file is faulty\n  (found only %d pixels instead of %d).\n",
+	   icount/pointsQE, ct_NPixels );
+  }
+
   // close file
 
   qefile.close();
+
+  // test QE
+
+  for(icount=0; icount< ct_NPixels; icount++){
+    for(i=0; i<pointsQE; i++){
+      if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
+	  QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
+	error( "read_pixels", "The quantum efficiency file is faulty\n  pixel %d, point %d  is % f, %f\n",
+	       icount, i, QE[icount][0][i], QE[icount][1][i] );
+      }
+    }
+  }
 
   // end
@@ -2875,16 +2133,4 @@
   } /* end if */
 
-  /* generate the ij coordinates */
-
-  for(i=0; i<pcam->inumpixels; i++){
-    pcam->dj[i] = pcam->dyc[i]/SIN60/dpsize;
-    pcam->di[i] = pcam->dxc[i]/dpsize - pcam->dj[i]*COS60;
-
-    //    fprintf(stdout, "%d %f %f %f %f %f\n", 
-    //	    i+1, pcam->di[i], pcam->dj[i], pcam->dxc[i], pcam->dyc[i],
-    //	    pcam->dpixsizefactor[i]);
-
-  }
-
   return(pcam->inumpixels);
 
@@ -2968,4 +2214,560 @@
           );
 }
+
+//------------------------------------------------------------
+// @name check_reflector_file                          
+//                                     
+// @desc check if a given reflector file has the right signature
+// @desc return TRUE or FALSE
+//
+// @date Mon Feb 14 16:44:21 CET 2000
+// @function @code 
+//------------------------------------------------------------
+
+int check_reflector_file(FILE *infile){
+
+  char Signature[20];           // auxiliary variable
+  char sign[20];                // auxiliary variable
+
+  strcpy(Signature, REFL_SIGNATURE);
+  
+  strcpy(sign, Signature);
+  
+  fread( (char *)sign, strlen(Signature), 1, infile);
+
+  if (strcmp(sign, Signature) != 0) {
+    cout << "ERROR: Signature of .rfl file is not correct\n";
+    cout << '"' << sign << '"' << '\n';
+    cout << "should be: " << Signature << '\n';
+    return(FALSE);
+  }
+
+  fread( (char *)sign, 1, 1, infile);
+
+  return(TRUE);
+
+}
+
+//------------------------------------------------------------
+// @name lin_interpol                          
+//                                     
+// @desc interpolate linearly between two points returning the 
+// @desc y-value of the result
+//
+// @date Thu Feb 17 11:31:32 CET 2000
+// @function @code 
+//------------------------------------------------------------
+
+float lin_interpol( float x1, float y1, float x2, float y2, float x){
+
+  if( (x2 - x1)<=0. ){ // avoid division by zero, return average
+    cout << "Warning: lin_interpol was asked to interpolate between two points with x1>=x2.\n";
+    return((y1+y2)/2.);
+  }
+  else{ // note: the check whether x1 < x < x2 is omitted for speed reasons
+    return((float) (y1 + (y2-y1)/(x2-x1)*(x-x1)) );
+  }
+}
+
+
+//------------------------------------------------------------
+// @name Photoelectron                          
+//                                     
+// @desc constructor for class Photoelectron
+//
+// @date Mon Feb 15 16:44:21 CET 2000
+// @function @code 
+//------------------------------------------------------------
+
+Photoelectron::Photoelectron(void){
+  iarrtime_ns = NOTIME;
+  ipixnum = -1;
+}
+
+//------------------------------------------------------------
+// @name produce_phes                          
+//                                     
+// @desc read the photons of an event, pixelize them and simulate QE
+// @desc return various statistics and the array of Photoelectrons
+//
+// @date Mon Feb 14 16:44:21 CET 2000
+// @function @code 
+//------------------------------------------------------------
+
+int produce_phes( FILE *sp, // the input file
+		  struct camera *cam, // the camera layout
+		  float minwl_nm, // the minimum accepted wavelength
+		  float maxwl_nm, // the maximum accepted wavelength
+		  class Photoelectron phe[iMAXNUMPHE], // the generated phes
+		  int *itotnphe, // total number of produced photoelectrons
+		  float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
+		  int *incph,    // total number of cph read
+		  float *tmin_ns,   // minimum arrival time of all phes
+		  float *tmax_ns    // maximum arrival time of all phes
+		  ){
+
+  static int i, k, ipixnum;
+  static float cx, cy, wl, qe, t;
+  static MCCphoton photon;
+  static float **qept;
+  static char flag[SIZE_OF_FLAGS + 1];
+  static float radius;
+
+  // reset variables
+
+  for ( i=0; i<cam->inumpixels; ++i ){
+    
+    nphe[i] = 0.0;
+    
+  }
+
+  *itotnphe = 0;
+  *incph = 0;
+  *tmin_ns = NOTIME; // very big 
+  *tmax_ns = -NOTIME; // very small
+
+  radius = cam->dxc[cam->inumpixels-1]
+    + 1.5*cam->dpixdiameter_cm*cam->dpixsizefactor[cam->inumpixels-1];
+
+  //- - - - - - - - - - - - - - - - - - - - - - - - - 
+  // read photons and "map" them into the pixels
+  //--------------------------------------------------      
+  
+  // initialize CPhoton
+  
+  photon.fill(0., 0., 0., 0., 0., 0., 0., 0.);
+  
+  // read the photons data
+  
+  fread ( flag, SIZE_OF_FLAGS, 1, sp );
+  
+  // loop over the photons
+  
+  while ( !isA( flag, FLAG_END_OF_EVENT ) ) {
+	  
+    memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
+    
+    fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
+	  
+    // increase number of photons
+	  
+    (*incph)++;
+	  	  
+    //
+    // Pixelization
+    //
+	    
+    cx = photon.get_x();
+    cy = photon.get_y(); 
+	  
+    // get wavelength
+	  
+    wl = photon.get_wl();
+
+    // cout << "wl " << wl << " radius " << sqrt(cx*cx + cy*cy) << "\n";
+	  
+    // check if photon has valid wavelength and is inside outermost camera radius
+	  
+    if( (wl > maxwl_nm) || (wl < minwl_nm) || (sqrt(cx*cx + cy*cy) > radius ) ){ 
+	    
+      // cout << " lost first check\n";
+
+      // read next CPhoton
+
+      fread ( flag, SIZE_OF_FLAGS, 1, sp );
+	    
+      // go to beginning of loop, the photon is lost
+      continue;
+	    
+    }
+
+    ipixnum = -1;
+
+    for(i=0; i<cam->inumpixels; i++){
+      if( bpoint_is_in_pix( cx, cy, i, cam) ){
+	ipixnum = i;
+	break;
+      }
+    }
+	   
+    if(ipixnum==-1){// the photon is in none of the pixels
+
+      // cout << " lost pixlization\n";
+
+      // read next CPhoton
+
+      fread ( flag, SIZE_OF_FLAGS, 1, sp );
+      
+      // go to beginning of loop, the photon is lost
+      continue;
+    }
+	  	  
+    //+++
+    // QE simulation
+    //---
+    
+    // set pointer to the QE table of the relevant pixel
+    	  
+    qept = (float **)QE[ipixnum];
+	  
+    // check if wl is inside table; outside the table, QE is assumed to be zero
+
+    if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
+
+      // cout << " lost\n"; 
+      
+      // read next Photon
+
+      fread ( flag, SIZE_OF_FLAGS, 1, sp );
+	    
+      // go to beginning of loop
+      continue;
+	    
+    }
+
+    // find data point in the QE table (-> k)
+
+    k = 1; // start at 1 because the condition was already tested for 0
+    while (k < pointsQE-1 && qept[0][k] < wl){
+      k++;
+    }
+
+    // calculate the qe value between 0. and 1.
+
+    qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
+    
+    // if random > quantum efficiency, reject it
+	        
+    if ( RandomNumber > qe ) {
+
+      // cout << " lost\n"; 
+	    
+      // read next Photon
+
+      fread ( flag, SIZE_OF_FLAGS, 1, sp );
+	    
+      // go to beginning of loop
+      continue;
+	    
+    }
+	  
+    //+++
+    // The photon has produced a photo electron
+    //---
+
+    // cout << " accepted\n";
+	  
+    // increment the number of photoelectrons in the relevant pixel
+	  
+    nphe[ipixnum] += 1.0;
+
+    t = photon.get_t() ; 
+
+    //    cout << " t " << t;
+	  
+    // find minimum and maximum arrival time
+
+    if(t < *tmin_ns){ 
+      *tmin_ns = t; // memorize time
+    }
+    if(t > *tmax_ns){ 
+      *tmax_ns = t; // memorize time
+    }
+
+    // store the new photoelectron
+
+    if(*itotnphe >= iMAXNUMPHE){
+      cout << "Error: Memory overflow. Event produces more than maximum\n";
+      cout << "       allowed number of photoelectrons (" << iMAXNUMPHE << ").\n";
+      return(1);
+    }
+
+    phe[*itotnphe].iarrtime_ns = (int)t;
+    phe[*itotnphe].ipixnum = ipixnum;
+	  
+    *itotnphe += 1;
+
+    // read next Photon
+	  
+    fread( flag, SIZE_OF_FLAGS, 1, sp );
+	  
+  }  // end while, i.e. found end of event
+	  
+  return(0);
+
+}
+
+
+//------------------------------------------------------------
+// @name produce_nsbrates
+//                                     
+// @desc read the starfield file, call produce_phes on it in,
+// @desc different wavebands, calculate the nsbrates
+//
+// @date Mon Feb 14 16:44:21 CET 2000
+// @function @code 
+//------------------------------------------------------------
+
+int produce_nsbrates( char *iname, // the starfield input file name
+		      struct camera *cam, // camera layout
+		      class Photoelectron phe[iMAXNUMPHE], // array for photoelectrons
+		      float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
+		                               // the NSB rates in phe/ns for each pixel 
+		      ){
+
+  int i, j, k, ii; // counters
+		      
+  static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
+					    WAVEBANDBOUND2,
+					    WAVEBANDBOUND3,
+					    WAVEBANDBOUND4,
+					    WAVEBANDBOUND5, 
+					    WAVEBANDBOUND6 };
+
+  FILE *infile;    // the input file
+  fpos_t fileposition; // marker on the input file
+  static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
+  static MCEventHeader evth; // the event header
+  static float nphe[iMAXNUMPIX];    // the number of photoelectrons in each pixel
+  int itnphe;   // total number of produced photoelectrons
+  int itotnphe;   // total number of produced photoelectrons after averaging
+  int incph;      // total number of cph read
+  float tmin_ns;  // minimum arrival time of all phes
+  float tmax_ns;  // maximum arrival time of all phes
+  float integtime_ns; // integration time from the starfield generator
+
+  // open input file 
+  
+  log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
+
+  infile = fopen( iname, "r" );
+  if ( infile == NULL ) 
+    error( SIGNATURE, "Cannot open starfield input file: %s\n", iname );
+  
+  // get signature, and check it
+  
+  if(check_reflector_file(infile)==FALSE){
+    exit(1);
+  }
+
+  // initialize flag
+  
+  strcpy( cflag, "                                        \0" );
+
+  // get flag
+    
+  fread( cflag, SIZE_OF_FLAGS, 1, infile );
+
+  if ( ! feof(infile)){ 
+
+    // reading .rfl file 
+
+    if(!isA( cflag, FLAG_START_OF_RUN )){
+      error( SIGNATURE, "Expected start of run flag, but found: %s\n", cflag );
+    }
+    else { // found start of run
+      
+      fread( cflag, SIZE_OF_FLAGS, 1, infile );
+      
+      if( isA( cflag, FLAG_START_OF_EVENT   )){ // there is a event
+	
+	// get MCEventHeader
+	
+	fread( (char*)&evth, evth.mysize(), 1, infile );
+	
+	integtime_ns = evth.get_energy();
+
+	// memorize where we are in the file
+
+	if (fgetpos( infile, &fileposition ) != 0){
+	  error( SIGNATURE, "Cannot position in file ...\n");
+	} 
+
+	// loop over the wavebands
+
+	for(i=0; i<iNUMWAVEBANDS; i++){
+
+	  // initialize the rate array
+
+	  for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
+	    rate_phepns[j][i] = 0.;
+	  }
+
+	  itotnphe = 0;
+
+	  // read the photons and produce the photoelectrons
+	  // - in order to average over the QE simulation, call the 
+          // production function iNUMNSBPRODCALLS times 
+
+	  for(ii=0; ii<iNUMNSBPRODCALLS; ii++){
+
+	    // position the file pointer to the beginning of the photons
+	    
+	    fsetpos( infile, &fileposition);
+
+	    // produce photoelectrons
+
+	    k = produce_phes( infile,
+			      cam,
+			      wl_nm[i],
+			      wl_nm[i+1],
+			      phe, // this is a dummy here
+			      &itnphe,
+			      nphe, // we want this!
+			      &incph,
+			      &tmin_ns,
+			      &tmax_ns );
+
+	    if( k != 0 ){ // non-zero returnvalue means error
+	      cout << "Exiting.\n";
+	      exit(1);
+	    }
+	  
+	    for(j = 0; j<cam->inumpixels; j++){ // loop over pixels
+	      rate_phepns[j][i] += nphe[j]/integtime_ns/(float)iNUMNSBPRODCALLS;
+	    }
+
+	    itotnphe += itnphe;
+
+	  } // end for(ii=0 ...
+
+	  fprintf(stdout, "Starfield, %6f - %6f nm: %d photoelectrons for %6f ns integration time\n",
+		  wl_nm[i], wl_nm[i+1], itotnphe/iNUMNSBPRODCALLS, integtime_ns);
+
+	} // end for(i=0 ...
+
+      }
+      else{
+	cout << "Starfield file contains no event.\nExiting.\n";
+	exit(1);
+      } // end if( isA ... event
+    } // end if ( isA ... run
+  }
+  else{
+    cout << "Starfield file contains no run.\nExiting.\n";
+    exit(1);
+  }
+
+  fclose( infile );
+
+  return(0);
+
+}
+
+
+//------------------------------------------------------------
+// @name produce_nsb_phes
+//                                     
+// @desc produce the photoelectrons from the nsbrates
+//
+// @date Thu Feb 17 17:10:40 CET 2000
+// @function @code 
+//------------------------------------------------------------
+
+int produce_nsb_phes( float *atmin_ns,
+		      float *atmax_ns,
+		      float theta_rad,
+		      struct camera *cam,
+		      float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
+		      float difnsbr_phepns[iMAXNUMPIX],
+		      float extinction[iNUMWAVEBANDS],
+		      float fnpx[iMAXNUMPIX],
+		      Photoelectron photo[iMAXNUMPHE],
+		      int *inphe,
+		      float base_mv[iMAXNUMPIX]){
+
+  float simtime_ns; // NSB simulation time
+  int i, j, k, ii; 
+  float zenfactor;  //  correction factor calculated from the extinction
+  int inumnsbphe;   //  number of photoelectrons caused by NSB
+
+  ii = *inphe; // avoid dereferencing
+		
+  // check if the arrival times are set; if not generate them
+
+  if(*atmin_ns == NOTIME){
+
+    *atmin_ns = 0.;
+    *atmax_ns = simtime_ns = SLICES*WIDTH_TIMESLICE;
+
+  }
+  else{ // extend the event time window by the given offsets
+
+    *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS;
+    *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS;
+
+    simtime_ns = *atmax_ns - *atmin_ns;
+
+    // make sure the simulated time is long enough for the FADC simulation
+
+    if(simtime_ns< SLICES*WIDTH_TIMESLICE){
+      *atmax_ns = *atmin_ns + SLICES*WIDTH_TIMESLICE;
+      simtime_ns = SLICES*WIDTH_TIMESLICE;
+    }
+      
+  }
+
+  // initialize baselines
+
+  for(i=0; i<cam->inumpixels; i++){
+    base_mv[i] = 0.;
+  }
+
+  // calculate baselines and generate phes
+
+  for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
+    
+    // calculate the effect of the atmospheric extinction 
+    
+    zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
+    
+    for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
+      
+      inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)  
+			  * simtime_ns );
+
+      base_mv[j] += inumnsbphe;
+
+      // randomize
+      
+      inumnsbphe =  ignpoi( inumnsbphe );
+      
+      // create the photoelectrons
+      
+      for(k=0; k < inumnsbphe; k++){
+
+	if(ii >= iMAXNUMPHE){
+	  cout << "Error: Memory overflow. NSB simulation produces more than maximum\n";
+	  cout << "       allowed number of photoelectrons (" << iMAXNUMPHE << ").\n";
+	  return(1);
+	}
+
+	photo[ii].iarrtime_ns = (int)(RandomNumber * simtime_ns + *atmin_ns );
+	photo[ii].ipixnum = j;
+	
+	// cout << "Created phe " << photo[ii].iarrtime_ns << " " 
+	//     << photo[ii].ipixnum << "\n";
+	
+	ii++; // increment total number of photoelectons
+	
+	fnpx[j] += 1.; // increment number of photoelectrons in this pixel
+
+      }
+      
+    } // end for(j=0 ...
+  } // end for(i=0 ...
+
+  // finish baseline calculation
+
+  for(i=0; i<cam->inumpixels; i++){
+    base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns;
+  }
+
+  *inphe = ii; // update the pointer
+
+  return(0);
+}
+
+
 // @endcode
 
@@ -2977,4 +2779,8 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.4  2000/01/25 08:36:23  petry
+// The pixelization in previous versions was buggy.
+// This is the first version with a correct pixelization.
+//
 // Revision 1.3  2000/01/20 18:22:17  petry
 // Found little bug which makes camera crash if it finds a photon
@@ -2996,218 +2802,4 @@
 // The "rootification" was done by Dirk Petry and Harald Kornmayer. 
 //
-// In the following you can see the README file of that version:
-//
-// ==================================================
-//
-// Fri Oct 22  1999   D.P.
-//
-// The MAGIC Monte Carlo System
-//
-// Camera Simulation Programme
-// ---------------------------
-//
-// 1) Description
-//
-// This version is the result of the fusion of H.K.'s
-// root_camera which is described below (section 2)
-// and another version by D.P. which had a few additional
-// useful features.
-//
-// The version compiles under Linux with ROOT 2.22 installed
-// (variable ROOTSYS has to be set).
-//
-// Compile as before simply using "make" in the root_camera
-// directory.
-//
-// All features of H.K.'s root_camera were retained.
-//
-// Additional features of this version are:
-//
-//   a) HBOOK is no longer used and all references are removed.
-//
-//   b) Instead of HBOOK, the user is given now the possibility of 
-//      having Diagnostic data in ROOT format as a complement
-//      to the ROOT Raw data.
-//
-//      This data is written to the file which is determined by
-//      the new input parameter "diag_file" in the camera parameter
-//      file.
-//
-//      All source code file belonging to this part have filenames
-//      starting with "MDiag".
-//
-//      The user can read the output file using the following commands
-//      in an interactive ROOT session:
-//
-//        	root [0] .L MDiag.so
-// 	root [1] new TFile("diag.root");
-// 	root [2] new TTreeViewer("T");
-// 	
-//      This brings up a viewer from which all variables of the
-//      TTree can be accessed and histogrammed. This example
-//      assumes that you have named the file "diag.root", that
-//      you are using ROOT version 2.22 or later and that you have
-//      the shared object library "MDiag.so" which is produced
-//      by the Makefile along with the executable "camera".
-//        
-//  !   The contents of the so-called diag file is not yet fixed.
-//  !   At the moment it is what J.C.G. used to put into the HBOOK
-//  !   ntuple. In future versions the moments calculation can be
-//  !   removed and the parameter list be modified correspondingly.
-//
-//   c) Now concatenated reflector files can be read. This is useful
-//      if you have run the reflector with different parameters but
-//      you want to continue the analysis with all reflector data
-//      going into ONE ROOT outputfile.
-//
-//      The previous camera version contained a bug which made reading 
-//      of two or more concatenated reflector files impossible.
-//
-//   d) The reflector output format was changed. It is now version
-//      0.4 .
-//      The change solely consists in a shortening of the flag
-//      definition in the file 
-//
-//            include-MC/MCCphoton.hxx  
-//
-//  !   IF YOU WANT TO READ REFLECTOR FORMAT 0.3, you can easily
-//  !   do so by recompiling camera with the previous version of
-//  !   include-MC/MCCphoton.hxx.
-//
-//      The change was necessary for saving space and better
-//      debugging. From now on, this format can be frozen.
-//
-//  !   For producing reflector output in the new format, you
-//  !   of course have to recompile your reflector with the
-//  !   new include-MC/MCCphoton.hxx .
-//
-//   e) A first version of the pixelization with the larger
-//      outer pixels is implemented. THIS IS NOT YET FULLY
-//      TESTED, but first rough tests show that it works
-//      at least to a good approximation.
-//
-//      The present version implements the camera outline
-//      with 18 "gap-pixels" and 595 pixels in total as
-//      shown in 
-//
-//         http://sarastro.ifae.es/internal/home/hardware/camera/numbering.ps
-//
-//      This change involved 
-//
-// 	(i) The file pixels.dat is no longer needed. Instead
-//           the coordinates are generated by the program itself
-//           (takes maybe 1 second). In the file 
-//
-// 		pixel-coords.txt
-//
-// 	  in the same directory as this README, you find a list
-//           of the coordinates generated by this new routine. It
-//           has the format
-//
-//               number   i   j   x  y  size-factor
-//
-//           where i and j are J.C.G.'s so called biaxis hexagonal
-//           coordinates (for internal use) and x and y are the
-//           coordinates of the pixel centers in the standard camera
-//           coordinate system in units of centimeters. The value
-//           of "size-factor" determines the linear size of the pixel
-//           relative to the central pixels. 
-//
-//         (ii) The magic.def file has two additional parameters
-//           which give the number of central pixels and the
-//           number of gap pixels
-//
-//         (iii) In camera.h and camera.cxx several changes were 
-//           necessary, among them the introduction of several
-//           new functions 
-//
-//      The newly suggested outline with asymmetric Winston cones
-//      will be implemented in a later version.
-//
-//   f) phe files can no longer be read since this contradicts
-//      our philosophy that the analysis should be done with other
-//      programs like e.g. EVITA and not with "camera" itself.
-//      This possibility was removed. 
-//
-//   g) ROOT is no longer invoked with an interactive interface.
-//      In this way, camera can better be run as a batch program and
-//      it uses less memory.
-//
-//   h) small changes concerning the variable "t_chan" were necessary in
-//      order to avoid segmentation faults: The variable is used as an
-//      index and it went sometimes outside the limits when camera
-//      was reading proton data. This is because the reflector files
-//      don't contain the photons in a chronological order and also
-//      the timespread can be considerably longer that the foreseen
-//      digitisation timespan. Please see the source code of camera.cxx
-//      round about line 1090.
-//
-//   j) several unused variables were removed, a few warning messages
-//      occur when you compile camera.cxx but these can be ignored at
-//      the moment.
-//
-// In general the program is of course not finished. It still needs
-// debugging, proper trigger simulation, simulation of the asymmetric
-// version of the outer pixels, proper NSB simulation, adaption of
-// the diag "ntuple" contents to our need and others small improvements.
-//
-// In the directory rfl-files there is now a file in reflector format 0.4
-// containing a single event produced by the starfiled adder. It has
-// a duration of 30 ns and represents the region around the Crab Nebula.
-// Using the enclosed input parameter file, camera should process this
-// file without problems.
-//
-// 2) The README for the previous version of root_camera
-//
-// README for a preliminary version of the 
-// root_camera program. 
-//
-// root_camera is based on the program "camera"of Jose Carlos
-// Gonzalez. It was changed in the way that only the pixelisation 
-// and the distibution of the phe to the FADCs works in a 
-// first version. 
-//
-// Using the #undef command most possibilities of the orignal 
-// program are switched of. 
-//
-// The new parts are signed by 
-//
-// - ROOT or __ROOT__ 
-//   nearly all  important codelines for ROOT output are enclosed 
-//   in structures like 
-//   #ifdef __ROOT__ 
-//   
-//     code 
-//
-//   #endif __ROOT__ 
-//
-//   In same case the new lines are signed by a comment with the word 
-//   ROOT in it. 
-//
-//   For timing of the pulse some variable names are changed. 
-//   (t0, t1, t  -->  t_ini, t_fin, t_1st, t_chan,...) 
-//   Look also for this changes. 
-//
-//   For the new root-file is also a change in readparm-files
-//
-//
-// - __DETAIL_TRIGGER__
-//
-//   This is for the implementation of the current work on trigger 
-//   studies. Because the class MTrigger is not well documented it 
-//   isn´t a part of this tar file. Only a dummy File exists. 
-//
-//
-//
-// With all files in the archive, the root_camera program should run. 
-//
-// A reflector file is in the directory rfl-files
-//
-// ==================================================
-//
-// From now on, use CVS for development!!!!
-//
-//
-//
 // Revision 1.3  1999/10/22 15:01:28  petry
 // version sent to H.K. and N.M. on Fri Oct 22 1999
