//!/////////////////////////////////////////////////////////////////////
//
// camera                
//
// @file        camera.cxx
// @title       Camera simulation
// @subtitle    Code for the simulation of the camera phase
// @desc        Code for the simulation of the camera of CT1 and MAGIC
// @author      J C Gonzalez
// @email       gonzalez@mppmu.mpg.de
// @date        Thu May  7 16:24:22 1998
//
//----------------------------------------------------------------------
//
// Created: Thu May  7 16:24:22 1998
// Author:  Jose Carlos Gonzalez
// Purpose: Program for reflector simulation
// Notes:   See files README for details
//    
//----------------------------------------------------------------------
//
// $RCSfile: camera.cxx,v $
// $Revision: 1.1.1.1 $
// $Author: harald $ 
// $Date: 1999-11-05 11:59:31 $
//
////////////////////////////////////////////////////////////////////////
// @tableofcontents @coverpage

//=-----------------------------------------------------------
//!@section Source code of |camera.cxx|.

/*!@"

  In this section we show the (commented) code of the program for the
  read-out of the output files generated by the simulator of the
  reflector, |reflector 0.3|.

  @"*/

//=-----------------------------------------------------------
//!@subsection Includes and Global variables definition.

//!@{

// includes for ROOT
// BEWARE: the order matters!

#include "TROOT.h"

#include "TFile.h"
#include "TTree.h"
#include "TBranch.h"

#include "MDiag.h"

#include "MTrigger.hxx"

#include "MRawEvt.h"
#include "MMcEvt.h"

/*!@" 

  All the defines are located in the file |camera.h|.

  @"*/

#include "camera.h"
//!@}

/*!@"

  The following set of flags are used in time of compilation. They do
  not affect directly the behaviour of the program at run-time
  (though, of course, if you disconnected the option for
  implementation of the Trigger logic, you will not be able to use any
  trigger at all. The 'default' values mean default in the sense of
  what you got from the server when you obtained this program.

  @"*/

//!@{

// flag for debugging (default: OFF )
#define __DEBUG__
#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__

//!@}

//=-----------------------------------------------------------
//!@subsection Definition of global variables.

/*!@"

  Now we define some global variables with data about the telescope, 
  such as "focal distance",  number of pixels/mirrors, 
  "size of the camera", and so on.

  @"*/

/*!@"

  Depending on the telescope we are using (CT1 or MAGIC), the 
  information stored in the definition file is different.
  The variable |ct_Type| has the value 0 when we use
  CT1, and 1 when we use MAGIC.

  @"*/

//!@{
static int   ct_Type;         //@< Type of telescope: 0:CT1, 1:MAGIC
//!@}

/*!@"

  And this is the information about the whole telescope.

  @"*/

//!@{

// parameters of the CT (from the CT definition file) 

////@: Focal distances [cm]
//static float *ct_Focal;       

//@: Mean Focal distances [cm]
static float ct_Focal_mean;   

//@: STDev. Focal distances [cm]
static float ct_Focal_std;    

//@: Mean Point Spread function [cm]
static float ct_PSpread_mean; 

//@: STDev. Point Spread function [cm]
static float ct_PSpread_std;  

//@: STDev. Adjustmente deviation [cm]
static float ct_Adjustment_std; 

//@: Radius of the Black Spot in mirror [cm]
static float ct_BlackSpot_rad;

//@: Radius of one mirror [cm]
static float ct_RMirror;      

//@: Camera width [cm]
static float ct_CameraWidth;  

//@: Pixel width [cm]
static float ct_PixelWidth;   

//@: ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(60)
static float ct_PixelWidth_corner_2_corner; 

//@: ct_PixelWidth_corner_2_corner / 2
static float ct_PixelWidth_corner_2_corner_half; 

//@: Number of mirrors
static int ct_NMirrors = 0;   

//@: Number of pixels
static int ct_NPixels;        

//@: Number of pixels
static int ct_NCentralPixels;        

//@: Number of pixels
static int ct_NGapPixels;        

//@: ct_Apot = ct_PixelWidth / 2
static float ct_Apot;          

//@: ct_2Apot = 2 * ct_Apot = ct_PixelWidth 
static float ct_2Apot;         

//@: name of the CT definition file to use
static char ct_filename[256];  

//@: list of showers to be skipped
static int *Skip;

//@: number of showers to be skipped
static int nSkip=0;

//@: flag: TRUE: data come from STDIN; FALSE: from file
static int Data_From_STDIN = FALSE;

//@: flag: TRUE: write all images to output; FALSE: only triggered showers
static int Write_All_Images = FALSE;

//@: flag: TRUE: write all data to output; FALSE: only triggered showers
static int Write_All_Data = FALSE;

//@: flag: TRUE: selection on the energy
static int Select_Energy = TRUE;

//@: Lower edge of the selected energy range (in GeV)
static float Select_Energy_le = 0.0; 

//@: Upper edge of the selected energy range (in GeV)
static float Select_Energy_ue = 100000.0; 

//!@}

/*!@"

  The following double-pointer is a 2-dimensional table with information 
  about each pixel. The routine read_pixels will generate
  the information for filling it using igen_pixel_coordinates().

  @"*/

//!@{
// Pointer to a tables/Arrays with information about the pixels
// 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
static float **pixary;  

//@: indexes of pixels neighbours of a given one
static int **pixneig;   

//@: number of neighbours a pixel have
static int *npixneig;   

//@: contents of the pixels (ph.e.)
static float *fnpix;    

//@: 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 ;

 
//!@}

/*!@"

  The following double-pointer is a 2-dimensional table with the
  Quantum Efficiency @$QE@$ of each pixel in the camera, as a function
  of the wavelength @$\lambda@$. The routine |read_pixels()| will read
  also this information from the file |qe.dat|.

  @"*/

//!@{
// Pointer to a table with QE, number of datapoints, and wavelengths

//@: table of QE
static float ***QE;

//@: number of datapoints for the QE curve
static int pointsQE;

//@: table of QE
static float *QElambda;
//!@}

/*!@"

  The following double-pointer is a 2-dimensional table with information 
  about each mirror in the dish. The routine |read_ct_file()| will read
  this information from the CT definition file.

  @"*/

//!@{
// Pointer to a table with the following info.:

static float **ct_data;       

/*
 *  TYPE=0  (CT1)
 *      i   s   rho   theta   x   y   z   thetan  phin  xn   yn   zn
 * 
 *       i : number of the mirror
 *       s : arc length [cm]
 *     rho : polar rho of the position of the center of the mirror [cm]
 *   theta : polar angle of the position of the center of the mirror [cm]
 *       x : x coordinate of the center of the mirror [cm]
 *       y : y coordinate of the center of the mirror [cm]
 *       z : z coordinate of the center of the mirror [cm]
 *  thetan : polar theta angle of the direction where the mirror points to
 *    phin : polar phi angle of the direction where the mirror points to
 *      xn : xn coordinate of the normal vector in the center (normalized)
 *      yn : yn coordinate of the normal vector in the center (normalized)
 *      zn : zn coordinate of the normal vector in the center (normalized)
 * 
 *  TYPE=1  (MAGIC)
 *      i  f   sx   sy   x   y   z   thetan  phin 
 * 
 *       i : number of the mirror
 *       f : focal distance of that mirror
 *      sx : curvilinear coordinate of mirror's center in X[cm]
 *      sy : curvilinear coordinate of mirror's center in X[cm]
 *       x : x coordinate of the center of the mirror [cm]
 *       y : y coordinate of the center of the mirror [cm]
 *       z : z coordinate of the center of the mirror [cm]
 *  thetan : polar theta angle of the direction where the mirror points to
 *    phin : polar phi angle of the direction where the mirror points to
 *      xn : xn coordinate of the normal vector in the center (normalized)
 *      yn : yn coordinate of the normal vector in the center (normalized)
 *      zn : zn coordinate of the normal vector in the center (normalized)
 */
//!@} 

/*!@"

  We define a table into where random numbers will be stored. 
  The routines used for random number generation are provided by
  |RANLIB| (taken from NETLIB, |www.netlib.org|), and by 
  the routine |double drand48(void)| (prototype defined in 
  |stdlib.h|) through the macro |RandomNumber| defined in 
  |camera.h|.

  @"*/

//!@{
// 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];  
//!@}

/*!@"

  The following are the set of parameters calculated for each image.
  The routines for their calculations are in |moments.cxx|.

  @"*/

//!@{
// parameters of the images

static Moments_Info *moments_ptr; 
static LenWid_Info *lenwid_ptr;

static float *maxs;
static int *nmaxs;
static float length, width, dist, xdist, azw, miss, alpha, *conc; 
static float phiasym, asymx, asymy;
static float charge, smax, maxtrigthr_phe;

//!@}

extern char FileName[];


//=-----------------------------------------------------------
// @subsection Main program.

//!@{

//++++++++++++++++++++++++++++++++++++++++
// MAIN PROGRAM 
//----------------------------------------

int main(int argc, char **argv) 
{

  //!@' @#### Definition of variables.
  //@'

  char inname[256];           //@< input file name
  char outname[256];          //@< output file name
  char datname[256];          //@< data (ASCII) output file name
  char diagname[256];         //@< diagnistic output file (ROOT format)
  char rootname[256] ;        //@< ROOT file name 

  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
  ofstream datafile;          //@< stream for the data file

  MCEventHeader mcevth;       //@< Event Header class (MC)
  MCCphoton cphoton;          //@< Cherenkov Photon class (MC)

  float thetaCT, phiCT;       //@< parameters of a given shower
  float thetashw, phishw;     //@< parameters of a given shower
  float coreD, coreX, coreY;  //@< core position and distance
  float impactD;              //@< impact parameter
  float l1, m1, n1;           //@< auxiliary variables
  float l2, m2, n2;           //@< auxiliary variables
  float num, den;             //@< auxiliary variables

  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 simulateNSB;            //@< Will we simulate NSB?
  float meanNSB;              //@< NSB mean value (per pixel)
  float qThreshold;           //@< Threshold value
  float qTailCut;             //@< Tail Cut value
  int nIslandsCut;            //@< Islands Cut value
  int countIslands;           //@< Will we count the islands?
  int anaPixels;
    
  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)
  float degTriggerZone;       //@< trigger area in the camera (radius, in deg.)

  float dtheta, dphi;         //@< deviations of CT from shower axis

  int still_in_loop = FALSE;

  char    Signature[20];

  float *image_data;
  int nvar, hidt;

  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()|.
  //@'

  int ch, errflg = 0;          //@< used by getopt

  /*!@'

    @#### Beginning of the program.
    
    We start with the main program. First we (could) make some
    presentation, and follows the reading of the parameters file (now
    from the |stdin|), the reading of the CT parameters file, and the
    creation of the output file, where the processed data will be
    stored.

  */

  //++
  // START
  //--

  // make unbuffered output

  cout.setf ( ios::stdio );

  // parse command line options (see reflector.h)
  
  parname[0] = '\0';

  optarg = NULL;
  while ( !errflg && ((ch = getopt(argc, argv, COMMAND_LINE_OPTIONS)) != -1) )
    switch (ch) {
    case 'f':
      strcpy(parname, optarg);
      break;
    case 'h':
      usage();
      break;
    default :
      errflg++;
    }
  
  // show help if error
  
  if ( errflg>0 )
    usage();

  // make some sort of presentation

  present();
  
  // read parameters file

  if ( strlen(parname) < 1 )
    readparam(NULL);
  else
    readparam(parname);

  // read data from file or from STDIN?

  Data_From_STDIN = get_data_from_stdin();

  // write all images, even those without trigger?

  Write_All_Images = get_write_all_images();

  // write all data (i.e., ph.e.s in pixels)

  Write_All_Data = get_write_all_data();

  // get filenames

  strcpy( inname, get_input_filename() );
  strcpy( outname, get_output_filename() );
  strcpy( datname, get_data_filename() );
  strcpy( diagname, get_diag_filename() );
  strcpy( rootname, get_root_filename() );
  strcpy( ct_filename, get_ct_filename() );

  // get different parameters of the simulation

  qThreshold = get_threshold();
  qTailCut = get_tail_cut();
  simulateNSB = get_nsb( &meanNSB );
  countIslands = get_islands_cut( &nIslandsCut );

  // get selections on the parameters
  
  Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
  
  // log filenames information

  log(SIGNATURE,
      "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n",
      "Filenames",
      "In", inname, 
      "Out", outname, 
      "Data", datname,
      "Diag", diagname,
      "ROOT",  rootname, 
      "CT", ct_filename);

  
  // log flags information

  log(SIGNATURE,
      "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
      "Flags",
      "Data_From_STDIN",   ONoff(Data_From_STDIN),  
      "Write_All_Images",  ONoff(Write_All_Images), 
      "Write_All_Data",    ONoff(Write_All_Data));
      
  // log parameters information
  
  log(SIGNATURE,
      "%s:\n\t%20s: %f\n\t%20s: %f\n\t%20s: %f %s\n\t%20s: %f %s\n",
      "Parameters",
      "q0 (Threshold)", qThreshold,
      "t0 (Tail-cut)", qTailCut,
      "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB),
      "i0 (Islands-cut)", nIslandsCut, ONoff(countIslands));
  
  // log selections
  
  log(SIGNATURE,
      "%s:\n\t%20s: %s  (%f:%f)\n",
      "Selections:",
      "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
  
  // set all random numbers seeds

  setall( get_seeds(0), get_seeds(1) );

  // get list of showers to evt. skip

  nSkip = get_nskip_showers();

  if (nSkip > 0) {
    Skip = new int[ nSkip ]; 
    get_skip_showers( Skip );

    log(SIGNATURE, "There are some showers to skip:\n");
    for (i=0; i<nSkip; ++i)
      log(SIGNATURE, "\tshower # %d\n", Skip[i]);
  }
  
  // read parameters from the ct.def file

  read_ct_file();

  // read pixels data

  read_pixels(&cam);

  // initialise ROOT

  TROOT simple("simple", "MAGIC Telescope Monte Carlo");

  // prepare ROOT tree for the diagnostic data 
  
  TFile *hfile;

  hfile = new TFile( diagname,"RECREATE", "MAGIC Telescope MC diagnostic data");
  
  
  // Create the ROOT Tree for the diagnostic data
  
  TTree *tree = new TTree("T","MAGIC Telescope MC diagnostic data");
  tree->SetAutoSave(100000000); 
  
  Int_t split = 1;
  Int_t bsize = 64000;
  MDiagEventobject  *event = 0;
  
  // Create one branch. If splitlevel is set, event is a superbranch
  // creating a sub branch for each data member of the Eventobject event.
  
  tree->Branch("event", "MDiagEventobject", &event, bsize, split);

#ifdef __ROOT__

  MRawEvt *Evt   = new MRawEvt() ; 
  MMcEvt  *McEvt = new MMcEvt (); 

  // initalize the ROOT file 
  //
  //     erzeuge ein Root file 
  //

  TFile outfile ( rootname , "RECREATE" ) ; 

  //
  //      create a Tree for the Event data stream 
  //

  TTree EvtTree("EvtTree","Events of Run");

  bsize=128000; split=1;

  EvtTree.Branch("MRawEvt","MRawEvt", 
  		 &Evt, bsize, split);

  EvtTree.Branch("MMcEvt","MMcEvt", 
  		 &McEvt, bsize, split);


  float  flli = 0. ; 
  unsigned short ulli = 0 ; 

#endif __ROOT__

  // for safety and for dimensioning image_data: count the elements in the 
  // diagnostic data branch

  i=0;
  i++; // "n"
  i++; // "primary"
  i++; // "energy"
  i++; // "cored"
  i++; // "impact"
  i++; // "xcore"
  i++; // "ycore"
  i++; // "theta"
  i++; // "phi"
  i++; // "deviations"
  i++; // "dtheta"
  i++; // "dphi"
  i++; // "trigger"
  i++; // "ncphs"
  i++; // "maxpassthr_phe"    
  i++; // "nphes"
  i++; // "nphes2"
  i++; // "length"
  i++; // "width"
  i++; // "dist"
  i++; // "xdist"
  i++; // "azw"
  i++; // "miss"
  i++; // "alpha"
  i++; // "conc2"
  i++; // "conc3"
  i++; // "conc4"
  i++; // "conc5"
  i++; // "conc6"
  i++; // "conc7"
  i++; // "conc8"
  i++; // "conc9"
  i++; // "conc10"
  i++; // "asymx"
  i++; // "asymy"
  i++; // "phiasym"

  nvar = i;
  image_data = new float[nvar];

  // set plate scale (deg/cm) and trigger area (deg)

  plateScale_cm2deg = ( ct_Type == 0 ) ? (0.244/2.1) : 0.030952381;

  if ( ! get_trigger_radius( &degTriggerZone ) )
    degTriggerZone = ( ct_Type == 0 ) ? (5.0) : (5.0);

  if ( ! get_correction( &fCorrection ) )
    fCorrection = 1.0;

  // number of pixels for parameters
    
  anaPixels = get_ana_pixels();
  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() ) 
      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';
    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 );
  datafile.open( datname );
  
  if ( outputfile.bad() ) 
    error( SIGNATURE, "Cannot open output file: %s\n", outname );
  
  // write signature

  outputfile.write( SIGNATURE, sizeof(SIGNATURE) );

  // initializes flag
  
  strcpy( flag, "                                        \0" );

  // allocate space for PMTs numbers of pixels
  
  fnpix = new float [ ct_NPixels ];
  fnpixclean = new float [ ct_NPixels ];

#ifdef __ROOT__
  
  fnslicesum = new float [ (2 * SLICES) ] ; 
  
  float slices   [ct_NPixels][ (2 * SLICES) ] ; 
  float slices2  [ct_NPixels][ SLICES ] ; 

  float trans    [ SLICES ] ; 
#endif __ROOT__ 

  
  moments_ptr = moments( anaPixels, NULL, NULL, 0.0, 1 );
        
  //!@' @#### Main loop.
  //@'

  //begin my version
                                           
  // get flag
    
  if ( Data_From_STDIN ) 
    cin.read( flag, SIZE_OF_FLAGS );
  else
    inputfile.read ( flag, SIZE_OF_FLAGS );

  // loop over the file

  still_in_loop = TRUE;

  while (
         ((! Data_From_STDIN) && (! inputfile.eof()))
         ||
         (Data_From_STDIN && still_in_loop)
         ) { 

    // reading .rfl files 
    if(!isA( flag, FLAG_START_OF_RUN )){
      error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag );
    }
    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 );

      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;
	log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
	
	// get MCEventHeader
	
	if ( Data_From_STDIN ) 
	  cin.read( (char*)&mcevth, mcevth.mysize() );
	else
	  mcevth.read( inputfile );
	
	// calculate core distance and impact parameter
	
	coreD = mcevth.get_core(&coreX, &coreY);
	
	// calculate impact parameter (shortest distance betwee the original
	// trajectory of the primary (assumed shower-axis) and the
	// direction where the telescope points to
	// 
	// we use the following equation, given that the shower core position
	// is (x1,y1,z1)=(x,y,0),the  trajectory is given by (l1,m1,n1),
	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
	// and (l2,m2,n2)
	//
	//               |                     |
	//               | x1-x2  y1-y2  z1-z2 |
	//               |                     |
	//             + |   l1     m1     n1  |
	//             - |                     |
	//               |   l2     m2     n2  |
	//               |                     |
	// dist = ------------------------------------        ( > 0 )
	//        [ |l1 m1|2   |m1 n1|2   |n1 l1|2 ]1/2
	//        [ |     |  + |     |  + |     |  ]
	//        [ |l2 m2|    |m2 n2|    |n2 l2|  ]
	//
	// playing a little bit, we get this reduced for in our case:
	//
	//
	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
	//         [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 - 
	//          2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
	
	// read the direction of the incoming shower
	
	thetashw = mcevth.get_theta();
	phishw = mcevth.get_phi();
	
	// calculate vector for shower
	
	l1 = sin(thetashw)*cos(phishw);
	m1 = sin(thetashw)*sin(phishw);
	n1 = cos(thetashw);
	
	// read the deviation of the telescope with respect to the shower
	
	mcevth.get_deviations ( &thetaCT, &phiCT );
	
	if ( (thetaCT == 0.) && (phiCT == 0.) ) {
	  
	  // CT was looking to the source (both lines are parallel)
	  // therefore, we calculate the impact parameter as the distance 
	  // between the CT axis and the core position
	  
	  impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
	  
	} else {
	  
	  // the shower comes off-axis
	  
	  // obtain with this the final direction of the CT
	  
	  thetaCT += thetashw;
	  phiCT   += phishw;
	  
	  // calculate vector for telescope
	  
	  l2 = sin(thetaCT)*cos(phiCT);
	  m2 = sin(thetaCT)*sin(phiCT);
	  n2 = cos(thetaCT);
	  
	  num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
	  den = (SQR(l1*m2 - l2*m1) + 
		 SQR(m1*n2 - m2*n1) + 
		 SQR(n1*l2 - n2*l1));
	  den = sqrt(den);
	  
	  impactD = fabs(num)/den;
	  
	  // fprintf(stderr, "[%f %f,%f %f] (%f %f %f) (%f %f %f) %f/%f = ",
	  //         thetashw, phishw, thetaCT, phiCT, l1, m1, n1, l2, m2, n2,
	  //         num, den);
	  
	}

	// 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__ 
	}

	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 (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(); 
	  
	  // get wavelength
	  
	  last_wl = wl;
	  wl = cphoton.get_wl();
	  
	  if ( wl < 1.0 )
	    break;
	  
	  if ( (wl > 600.0) || (wl < 290.0) )
	    break;
	  
	  // check if photon is inside outermost camera radius

	  if(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;
	  
	  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];
	  
	  }
	  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 ) {
	    
	    // 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;
	
	// skip it ?
	
	for ( i=0; i<nSkip; ++i ) {
	  if (Skip[i] == (nshow+ntshow)) {
	    i = -1;
	    break;
	  }
	}
	
	// if after the previous loop, the exit value of i is -1
	// then the shower number is in the list of showers to be
	// skipped
	
	if (i == -1) {
	  log(SIGNATURE, "\t\tskipped!\n");
	  continue;
	}
	
	/*!@'
	  
	  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
	
	if ( Select_Energy ) {
	  if (( mcevth.get_energy() < Select_Energy_le ) ||
	      ( mcevth.get_energy() > Select_Energy_ue )) {
	    log(SIGNATURE, "select_energy: shower rejected.\n");
	    continue;
	  }
	}
	
#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 ) ; 
	  
	}
      }
      
      //
      //   
      //
      
      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 __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__
	
	/*!@'
	  
	  @#### 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;
	
	if ( trigger == TRUE ) {
	  
	  itrigger = i;
	  ++ntrigger;
	  
	  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__
	  
	  //!@' @#### Islands algorithm.
	  //@'
	  
	  //++
	  // islands counting, and cleanning
	  //--
	  
	  if ( countIslands )
	    do_islands( ct_NPixels, fnpixclean, pixneig, npixneig, 
			countIslands, nIslandsCut);
	  
#endif // __ISLANDS__
	  
#ifdef __MOMENTS__
	  
	  //!@' @#### Calculation of parameters of the image.
	  //@'
	  
	  //++
	  // moments calculation
	  //--
	  
	  // calculate moments and other things

	  moments_ptr = moments( anaPixels, fnpixclean, pixary,
                               plateScale_cm2deg, 0 );
	  
	  charge = moments_ptr->charge ;
	  smax   = moments_ptr->smax   ;
	  maxs   = moments_ptr->maxs   ;
	  nmaxs  = moments_ptr->nmaxs  ;
	  length = moments_ptr->length ;
	  width  = moments_ptr->width  ;
	  dist   = moments_ptr->dist   ;
	  xdist  = moments_ptr->xdist  ;
	  azw    = moments_ptr->azw    ;
	  miss   = moments_ptr->miss   ;
	  alpha  = moments_ptr->alpha  ;
	  conc   = moments_ptr->conc   ;
	  asymx  = moments_ptr->asymx  ;
	  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
	  
	  event = new MDiagEventobject();
	  
	  i=0;
	  image_data[i] =	event->n       = hidt/10; i++;
	  image_data[i] =       event->primary = mcevth.get_primary(); i++;
	  image_data[i] =	event->energy  = mcevth.get_energy(); i++;
	  image_data[i] =	event->cored   = coreD;  i++;
	  image_data[i] =	event->impact  = impactD; i++;
	  image_data[i] =	event->xcore   = coreX; i++;
	  image_data[i] =	event->ycore   = coreY; i++;
	  image_data[i] =	event->theta   = mcevth.get_theta(); i++;
	  image_data[i] =	event->phi     = mcevth.get_phi(); i++;
	  image_data[i] =       event->deviations = mcevth.get_deviations (&dtheta, &dphi); i++;
	  image_data[i] =       event->dtheta  = dtheta; i++;
	  image_data[i] =       event->dphi    = dphi; i++;
	  image_data[i] =       event->trigger = trigger; i++;
	  image_data[i] =	event->ncphs   = ncph; i++;
	  image_data[i] =	event->maxpassthr_phe   = maxtrigthr_phe; i++;
	  image_data[i] =	event->nphes   = charge; i++;
	  image_data[i] =	event->nphes2  = smax; i++;
	  image_data[i] =	event->length  = length; i++;
	  image_data[i] =	event->width   = width; i++;
	  image_data[i] =	event->dist    = dist; i++;
	  image_data[i] =	event->xdist   = xdist; i++;
	  image_data[i] =	event->azw     = azw; i++;
	  image_data[i] =	event->miss    = miss; i++;
	  image_data[i] =	event->alpha   = alpha; i++;
	  image_data[i] =	event->conc2   = conc[0]; i++;
	  image_data[i] =	event->conc3   = conc[1]; i++;
	  image_data[i] =	event->conc4   = conc[2]; i++;
	  image_data[i] =	event->conc5   = conc[3]; i++;
	  image_data[i] =	event->conc6   = conc[4]; i++;
	  image_data[i] =	event->conc7   = conc[5]; i++;
	  image_data[i] =	event->conc8   = conc[6]; i++;
	  image_data[i] =	event->conc9   = conc[7]; i++;
	  image_data[i] =	event->conc10  = conc[8]; i++;
	  image_data[i] =	event->asymx   = asymx; i++;
	  image_data[i] =	event->asymy   = asymy; i++;
	  image_data[i] =	event->phiasym = phiasym; i++;
	  
	  // there should be "nvar" variables
	  
	  if ( i != nvar ) 
	    error( SIGNATURE, "Wrong entry number for diagnostic data.\n" );
	  
	  tree->Fill();
	  delete event;
	  
	  // put information in the data file, 
	  
	  datafile << ntrigger;
	  for(i=0;i<nvar;i++) {
	    datafile << ' ' << image_data[i];
	  }
	  
	  
#endif // __MOMENTS__
	  
	  
	  // revert the fnpixclean matrix into fnpix
	  // (now we do this, but maybe in a future we want to
	  // use both fnpix and fnpixclean for different things
	  
	  memcpy( fnpix, fnpixclean, sizeof(float) * ct_NPixels );
	  
	  // put this information in the data file, 
	  
	  if ( Write_All_Data ) {
	    datafile << ' ' << -9999;
	    for ( i=0; i<ct_NPixels; ++i )
	      datafile << ' ' << fnpix[i];
	  }
	  
	  datafile << endl;  
	  
	  mcevth.set_trigger( TRUE );
	  
	  log(SIGNATURE, "TRIGGER\n");
	  
	} else { // ( trigger == FALSE )
	  
	  event = new MDiagEventobject();
	  
	  i=0;
	  image_data[i] =	event->n       = hidt/10; i++;
	  image_data[i] =       event->primary = mcevth.get_primary(); i++;
	  image_data[i] =	event->energy  = mcevth.get_energy(); i++;
	  image_data[i] =	event->cored   = coreD = mcevth.get_core(&coreX, &coreY); i++;
	  image_data[i] =	event->impact  = coreD; i++;
	  image_data[i] =	event->xcore   = coreX; i++;
	  image_data[i] =	event->ycore   = coreY; i++;
	  image_data[i] =	event->theta   = mcevth.get_theta(); i++;
	  image_data[i] =	event->phi     = mcevth.get_phi(); i++;
	  image_data[i] =       event->deviations = mcevth.get_deviations(&dtheta, &dphi); i++;
	  image_data[i] =       event->dtheta = dtheta; i++;
	  image_data[i] =       event->dphi = dphi; i++;
	  image_data[i] =       event->trigger = trigger; i++;
	  image_data[i] =	event->ncphs   = ncph; i++;
	  image_data[i] =	event->maxpassthr_phe   = maxtrigthr_phe; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  image_data[i] =	-1.; i++;
	  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" );
	  
	  tree->Fill();
	  delete event;
	  
	  // put this information in the data file, 

	  if ( Write_All_Data ) {

	    datafile << ntrigger;
	    for ( i=0; i<nvar; ++i )
	      datafile << ' ' << image_data[i];
	    
	    datafile << -9999;
	    for ( i=0; i<ct_NPixels; ++i )
	      datafile << ' ' << fnpix[i];
	    
	    datafile << endl;  
	  }
	  
	  mcevth.set_trigger( FALSE );
	  
	} // trigger == FALSE 
        
#endif // __TRIGGER__

	//!@' @#### Save data.
	//@'
	
	//++++++++++++++++++++++++++++++++++++++++++++++++++
	// we now have all information we want
	// the only thing we must do now is writing it to 
	// the output file
	//--------------------------------------------------

	//++ 
	// 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 );
	
      } // 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 );
      }
      else { // found end of run
	ntshow += nshow;
	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 );
	
	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
	  log(SIGNATURE, "End of file . . .\n");
	  still_in_loop  = FALSE;
	  
	  if ((! Data_From_STDIN) && (! inputfile.eof())){
	    
	    // 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';
	      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 );
    } // 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
      //--

      EvtTree.Write() ; 
      outfile.Write() ;
      outfile.Close() ; 

#endif __ROOT__
              
  // close input file
  
  ntcph += ncph;
  log( SIGNATURE, "%d event(s), with a total of %d C.photons\n", 
       ntshow, ntcph );
  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();
  datafile.close();

  hfile->Write();

  hfile->Close();

#ifdef __DETAIL_TRIGGER__ 
  // Output of Trigger statistics
  //

  Trigger.PrintStat() ; 
#endif __DETAIL_TRIGGER__ 

  // program finished

  log( SIGNATURE, "Done.\n");
  
  return( 0 );
}
//!@}

// @T \newpage

//!@subsection Functions definition.

//!-----------------------------------------------------------
// @name present
//
// @desc Make some presentation
//
// @date Sat Jun 27 05:58:56 MET DST 1998
//------------------------------------------------------------
// @function

//!@{
void 
present(void)
{
  cout << "##################################################\n"
       <<  SIGNATURE << '\n' << '\n'
       << "Processor of the reflector output\n"
       << "J C Gonzalez, Jun 1998\n"
       << "##################################################\n\n"
       << flush ;
}
//!@}


//!-----------------------------------------------------------
// @name usage
//
// @desc show help
//
// @date Tue Dec 15 16:23:30 MET 1998
//------------------------------------------------------------
// @function 

//!@{
void 
usage(void)
{
  present();
  cout << "\nusage ::\n\n"
       << "\t camera "
       << " [ -@ paramfile ] "
       << " [ -h ] "
       << "\n\n or \n\n"
       << "\t camera < paramfile"
       << "\n\n";
  exit(0);
}
//!@}


//!-----------------------------------------------------------
// @name log                             
//                                   
// @desc function to send log information
//
// @var    funct  Name of the caller function
// @var    fmt    Format to be used (message)
// @var    ...    Other information to be shown
//
// @date Sat Jun 27 05:58:56 MET DST 1998
//------------------------------------------------------------
// @function  

//!@{
void
log(const char *funct, char *fmt, ...)
{
  va_list args;
  
  //  Display the name of the function that called error
  printf("[%s]: ", funct);
  
  // Display the remainder of the message
  va_start(args, fmt);
  vprintf(fmt, args);
  va_end(args);
}
//!@}


//!-----------------------------------------------------------
// @name error                                                    
//                                                           
// @desc function to send an error message, and abort the program
//
// @var    funct  Name of the caller function
// @var    fmt    Format to be used (message)
// @var    ...    Other information to be shown
//
// @date Sat Jun 27 05:58:56 MET DST 1998
//------------------------------------------------------------
// @function  

//!@{
void
error(const char *funct, char *fmt, ...)
{
  va_list args;

  //  Display the name of the function that called error
  fprintf(stderr, "ERROR in %s: ", funct);

  // Display the remainder of the message
  va_start(args, fmt);
  vfprintf(stderr, fmt, args);
  va_end(args);

  perror(funct);

  abort();
}
//!@}


//!-----------------------------------------------------------
// @name isA                               
//                                             
// @desc returns TRUE(FALSE), if the flag is(is not) the given
//
// @var    s1     String to be searched
// @var    flag   Flag to compare with string s1
// @return TRUE: both strings match; FALSE: oth.
//
// @date Wed Jul  8 15:25:39 MET DST 1998
//------------------------------------------------------------
// @function  

//!@{
int 
isA( char * s1, const char * flag ) {
  return ( (strncmp((char *)s1, flag, SIZE_OF_FLAGS)==0) ? 1 : 0 ); 
}
//!@}


//!-----------------------------------------------------------
// @name read_ct_file           
//                          
// @desc read CT definition file
//
// @date Sat Jun 27 05:58:56 MET DST 1998
//------------------------------------------------------------
// @function  

//!@{
void
read_ct_file(void)
{
  char line[LINE_MAX_LENGTH];    //@< line to get from the ctin
  char token[ITEM_MAX_LENGTH];   //@< a single token
  int i, j;                      //@< dummy counters

  log( "read_ct_file", "start.\n" );

  ifstream ctin ( ct_filename );

  if ( ctin.bad() ) 
    error( "read_ct_file", 
           "Cannot open CT def. file: %s\n", ct_filename );
  
  // loop till the "end" directive is reached

  while (!ctin.eof()) {          

    // get line from stdin

    ctin.getline(line, LINE_MAX_LENGTH);

    // look for each item at the beginning of the line

    for (i=0; i<=define_mirrors; i++) 
      if (strstr(line, CT_ITEM_NAMES[i]) == line)
        break;
    
    // if it is not a valid line, just ignore it

    if (i == define_mirrors+1) 
      continue;
    
    // case block for each directive

    switch ( i ) {

    case type:                // <type of telescope> (0:CT1 ¦ 1:MAGIC)
      
      // get focal distance

      sscanf(line, "%s %d", token, &ct_Type);

      log( "read_ct_file", "<Type of Telescope>: %s\n", 
           ((ct_Type==0) ? "CT1" : "MAGIC") );

      break;

    case focal_distance:      // <focal distance> [cm]
      
      // get focal distance

      sscanf(line, "%s %f", token, &ct_Focal_mean);

      log( "read_ct_file", "<Focal distance>: %f cm\n", ct_Focal_mean );

      break;

    case focal_std:           // s(focal distance) [cm]
      
      // get focal distance

      sscanf(line, "%s %f", token, &ct_Focal_std);

      log( "read_ct_file", "s(Focal distance): %f cm\n", ct_Focal_std );

      break;

    case point_spread:        // <point spread> [cm]
      
      // get point spread

      sscanf(line, "%s %f", token, &ct_PSpread_mean);

      log( "read_ct_file", "<Point spread>: %f cm\n", ct_PSpread_mean );

      break;

    case point_std:           // s(point spread) [cm]
      
      // get point spread

      sscanf(line, "%s %f", token, &ct_PSpread_std);

      log( "read_ct_file", "s(Point spread): %f cm\n", ct_PSpread_std );

      break;

    case adjustment_dev:      // s(adjustment_dev) [cm]
      
      // get point spread

      sscanf(line, "%s %f", token, &ct_Adjustment_std);

      log( "read_ct_file", "s(Adjustment): %f cm\n", ct_Adjustment_std );

      break;

    case black_spot:          // radius of the black spot in the center [cm]
      
      // get black spot radius

      sscanf(line, "%s %f", token, &ct_BlackSpot_rad);

      log( "read_ct_file", "Radius of the black spots: %f cm\n", 
           ct_BlackSpot_rad);

      break;

    case r_mirror:            // radius of the mirrors [cm]
      
      // get radius of mirror

      sscanf(line, "%s %f", token, &ct_RMirror);

      log( "read_ct_file", "Radii of the mirrors: %f cm\n", ct_RMirror );

      break;

    case n_mirrors:           // number of mirrors
      
      // get the name of the output_file from the line

      sscanf(line, "%s %d", token, &ct_NMirrors);

      log( "read_ct_file", "Number of mirrors: %d\n", ct_NMirrors );

      break;

    case camera_width:        // camera width [cm]
      
      // get the name of the ct_file from the line

      sscanf(line, "%s %f", token, &ct_CameraWidth);

      log( "read_ct_file", "Camera width: %f cm\n", ct_CameraWidth );

      break;

    case n_pixels:           // number of pixels
      
      // get the name of the output_file from the line

      sscanf(line, "%s %d", token, &ct_NPixels);

      log( "read_ct_file", "Number of pixels: %d\n", ct_NPixels );

      break;

    case n_centralpixels:           // number of central pixels
      
      // get the name of the output_file from the line

      sscanf(line, "%s %d", token, &ct_NCentralPixels);

      log( "read_ct_file", "Number of central pixels: %d\n", ct_NCentralPixels );

      break;

    case n_gappixels:           // number of gap pixels
      
      // get the name of the output_file from the line

      sscanf(line, "%s %d", token, &ct_NGapPixels);

      log( "read_ct_file", "Number of gap pixels: %d\n", ct_NGapPixels );

      break;

    case pixel_width:         // pixel width [cm]
      
      // get the name of the ct_file from the line

      sscanf(line, "%s %f", token, &ct_PixelWidth);

      ct_PixelWidth_corner_2_corner = ct_PixelWidth / cos(RAD(30.0));
      ct_PixelWidth_corner_2_corner_half =
        ct_PixelWidth_corner_2_corner * 0.50;
      ct_Apot = ct_PixelWidth / 2;
      ct_2Apot = ct_Apot * 2.0;

      log( "read_ct_file", "Pixel width: %f cm\n", ct_PixelWidth );

      break;

    case define_mirrors:      // read table with the parameters of the mirrors

      log( "read_ct_file", "Table of mirrors data:\n" );

      // check whether the number of mirrors was already set

      if ( ct_NMirrors == 0 )
        error( "read_ct_file", "NMirrors was not set.\n" );
      
      // allocate memory for paths list

      log( "read_ct_file", "Allocating memory for ct_data\n" );

      ct_data = new float*[ct_NMirrors];

      for (i=0; i<ct_NMirrors; i++) 
        ct_data[i] = new float[CT_NDATA];

      // read data

      log( "read_ct_file", "Reading mirrors data...\n" );

      for (i=0; i<ct_NMirrors; i++)
        for (j=0; j<CT_NDATA; j++)
          ctin >> ct_data[i][j];

      break;

    } // switch ( i ) 

  } // end while

  // end

  log( "read_ct_file", "done.\n" );

  return;
}
//!@}


//!-----------------------------------------------------------
// @name read_pixels  
//                          
// @desc read pixels data
//
// @date Fri Mar 12 16:33:34 MET 1999
//------------------------------------------------------------
// @function

//!@{
void 
read_pixels(struct camera *pcam)
{
  ifstream qefile;
  char line[LINE_MAX_LENGTH];
  int n, i, j, k;
  float qe;

  //------------------------------------------------------------
  // first, pixels' coordinates

  pcam->inumpixels = ct_NPixels;
  pcam->inumcentralpixels = ct_NCentralPixels;
  pcam->inumgappixels = ct_NGapPixels;
  pcam->inumbigpixels = ct_NPixels - ct_NCentralPixels - ct_NGapPixels;
  pcam->dpixdiameter_cm =  ct_PixelWidth; 

  // 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];
  for ( i=0; i<2*ct_NCentralPixels; ++i ) 
    pixary[i] = new float[2];

  pixneig = new int* [ct_NCentralPixels];
  for ( i=0; i<ct_NCentralPixels; ++i ) {
    pixneig[i] = new int[6];
    for ( j=0; j<6; ++j ) 
      pixneig[i][j] = -1;
  }

  npixneig = new int[ct_NCentralPixels];
  for ( i=0; i<ct_NCentralPixels; ++i ) 
    npixneig[i] = 0;

  // generate all coordinates

  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
  
#ifdef __DEBUG__
  for ( n=0 ; n<ct_NPixels ; ++n ) {
    cout << "Para el pixel " << n << ": ";	
    for ( i=n+1 ; (i<ct_NPixels)&&(npixneig[n]<6) ; ++i) {
      if ( pixels_are_neig(n,i) == TRUE ) {
        pixneig[n][npixneig[n]] = i;
        pixneig[i][npixneig[i]] = n;
        cout << i << ' ';
        ++npixneig[n];
        ++npixneig[i];
      }
    }
    cout << endl << flush;
  }
#else // ! __DEBUG__
  for ( n=0 ; n<ct_NCentralPixels ; ++n ) 
    for ( i=n+1 ; (i<ct_NCentralPixels)&&(npixneig[n]<6) ; ++i) 
      if ( pixels_are_neig(n,i) == TRUE ) {
        pixneig[n][npixneig[n]] = i;
        pixneig[i][npixneig[i]] = n;
        ++npixneig[n];
        ++npixneig[i];
      }
#endif // ! __DEBUG__
  
#ifdef __DEBUG__
  for ( n=0 ; n<ct_NPixels ; ++n ) {
    cout << n << ':';
    for ( j=0; j<npixneig[n]; ++j) 
      cout << ' ' << pixneig[n][j];
    cout << endl << flush;
  }
#endif // __DEBUG__  

  //------------------------------------------------------------
  // second, pixels' QE

  // try to open the file

  log("read_pixels", "Openning the file \"%s\" . . .\n", QE_FILE);
  
  qefile.open( QE_FILE );
  
  // if it is wrong or does not exist, exit
  
  if ( qefile.bad() )
    error( "read_pixels", "Cannot open \"%s\". Exiting.\n", QE_FILE );
  
  // read file

  log("read_pixels", "Reading data . . .\n");

  i=-1;

  while ( ! qefile.eof() ) {          

    // get line from the file

    qefile.getline(line, LINE_MAX_LENGTH);

    // skip if comment

    if ( *line == '#' )
      continue;

    // if it is the first valid value, it is the number of QE data points

    if ( i < 0 ) {

      // get the number of datapoints 

      sscanf(line, "%d", &pointsQE);
      
      // allocate memory for the table of QEs
      
      QE = new float ** [ct_NPixels];

      for ( i=0; i<ct_NPixels; ++i ) {
        QE[i] = new float * [2];
        QE[i][0] = new float[pointsQE];
        QE[i][1] = new float[pointsQE];
      }
      
      QElambda = new float [pointsQE];

      for ( i=0; i<pointsQE; ++i ) {
        qefile.getline(line, LINE_MAX_LENGTH);
        sscanf(line, "%f", &QElambda[i]);
      }

      i=0;

      continue;
    }

    // get the values (num-pixel, num-datapoint, QE-value)
    
    sscanf(line, "%d %d %f", &i, &j, &qe);

    if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
         ((j-1) < pointsQE)   && ((j-1) > -1) ) {
      QE[i-1][0][j-1] = QElambda[j-1];
      QE[i-1][1][j-1] = qe;
    }

  }

  // close file

  qefile.close();

  // end

  log("read_pixels", "Done.\n");

}
//!@}


//!-----------------------------------------------------------
// @name pixels_are_neig                        
//                                             
// @desc check whether two pixels are neighbours
//
// @var pix1      Number of the first pixel
// @var pix2      Number of the second pixel
// @return        TRUE: both pixels are neighbours; FALSE: oth.
//
// @date Wed Sep  9 17:58:37 MET DST 1998
//------------------------------------------------------------
// @function  

//!@{
int
pixels_are_neig(int pix1, int pix2)
{ 
  if ( sqrt(SQR( pixary[pix1][0] - pixary[pix2][0] ) +
            SQR( pixary[pix1][1] - pixary[pix2][1] ) ) 
       > ct_PixelWidth_corner_2_corner ) 
    return ( FALSE );
  else
    return ( TRUE );
}
//!@}

//!-----------------------------------------------------------
// @name igen_pixel_coordinates
//                                             
// @desc generate the pixel center coordinates
//
// @var *pcam     structure camera containing all the
//                camera information
// @return        total number of pixels
//
// DP
//
// @date Thu Oct 14 10:41:03 CEST 1999
//------------------------------------------------------------
// @function  

//!@{
/******** igen_pixel_coordinates() *********************************/

int igen_pixel_coordinates(struct camera *pcam) { 
            /* generate pixel coordinates, return value is number of pixels */

  int i, itot_inside_ring, iN, in, ipixno, iring_no, ipix_in_ring, isegment;
  float fsegment_fract;
  double dtsize;
  double dhsize;
  double dpsize;
  double dxfirst_pix;
  double dyfirst_pix;
  double ddxseg1, ddxseg2, ddxseg3, ddxseg4, ddxseg5, ddxseg6;
  double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5, ddyseg6;
  

  double dstartx, dstarty;   /* for the gap pixels and outer pixels */
  int j, nrow;

  dpsize = pcam->dpixdiameter_cm;
  dtsize = dpsize * sqrt(3.) / 2.;
  dhsize = dpsize / 2.;

  /* Loop over central pixels to generate co-ordinates  */

  for(ipixno=1; ipixno <= pcam->inumcentralpixels; ipixno++){

    /* Initialise variables. The central pixel = ipixno 1 in ring iring_no 0 */

    pcam->dpixsizefactor[ipixno] = 1.;

    in = 0;

    i = 0;
    itot_inside_ring = 0;
    iring_no = 0;

    /* Calculate the number of pixels out to and including the ring containing pixel number */
    /* ipixno e.g. for pixel number 17 in ring number 2 itot_inside_ring = 19 */

    while (itot_inside_ring == 0){
      
      iN = 3*(i*(i+1)) + 1;
      
      if (ipixno <= iN){
	iring_no = i;
	itot_inside_ring = iN;
      }
      
      i++;
    }
    
    
    /* Find the number of pixels which make up ring number iring_no e.g. ipix_in_ring = 6 for ring 1 */    
        
    ipix_in_ring = 0;
    for (i = 0; i < iring_no; ++i){

      ipix_in_ring = ipix_in_ring + 6;
    }

    /* The camera is viewed as 6 radial segments ("pie slices"). Knowing the number of pixels in its */
    /* ring calculate which segment the pixel ipixno is in. Then find how far across this segment it is */
    /* as a fraction of the number of pixels in this sixth of the ring (ask SMB). */
        
    isegment = 0;
    fsegment_fract = 0.;
    if (iring_no > 0) {
      
      isegment = (int)((ipixno - itot_inside_ring + ipix_in_ring - 0.5) / iring_no + 1); /* integer division ! numbering starts at 1 */
      
      fsegment_fract = (ipixno - (itot_inside_ring - ipix_in_ring)) - ((isegment-1)*iring_no) - 1 ;
      
    }

    /* the first pixel in each ring lies on the positive x axis at a distance dxfirst_pix = iring_no * the */
    /* pixel width (flat to flat) dpsize. */
        
    dxfirst_pix = dpsize*iring_no;
    dyfirst_pix = 0.;

    /* the vector between the first and last pixels in a segment n is (ddxsegn, ddysegn) */

    ddxseg1 = - dhsize*iring_no;
    ddyseg1 = dtsize*iring_no;
    ddxseg2 = -dpsize*iring_no;
    ddyseg2 = 0.;
    ddxseg3 = ddxseg1;
    ddyseg3 = -ddyseg1;
    ddxseg4 = -ddxseg1;
    ddyseg4 = -ddyseg1;
    ddxseg5 = -ddxseg2;
    ddyseg5 = 0.;
    ddxseg6 = -ddxseg1;
    ddyseg6 = ddyseg1;
    
    /* to find the position of pixel ipixno take the position of the first pixel in the ring and move */
    /* anti-clockwise around the ring by adding the segment to segment vectors. */

    switch (isegment) {
      
    case 0:

      pcam->dxc[ipixno-1] = 0.;
      pcam->dyc[ipixno-1] = 0.; 

    case 1: 
      pcam->dxc[ipixno-1] = dxfirst_pix - dhsize*fsegment_fract;
      pcam->dyc[ipixno-1] = dyfirst_pix + dtsize*fsegment_fract;
      
      break;
      
    case 2:
      
      pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 - dpsize*fsegment_fract;
      pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + 0.;
      
      break;
      
    case 3:
      
      pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 - dhsize*fsegment_fract;
      pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 - dtsize*fsegment_fract;
      
      break;
      
    case 4:
      
      pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + dhsize*fsegment_fract;
      pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 - dtsize*fsegment_fract;
      
      break;
      
    case 5:
      
      pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + dpsize*fsegment_fract;
      pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + 0.;
    
      break;
      
    case 6:
      
      pcam->dxc[ipixno-1] = dxfirst_pix + ddxseg1 + ddxseg2 + ddxseg3 + ddxseg4 + ddxseg5 + dhsize*fsegment_fract;
      pcam->dyc[ipixno-1] = dyfirst_pix + ddyseg1 + ddyseg2 + ddyseg3 + ddyseg4 + ddyseg5 + dtsize*fsegment_fract;
      
      break;
      
    default: 

      fprintf(stderr, "ERROR: problem in coordinate generation for pixel %d\n", ipixno);
      return(0);
      
    } /* end switch */

  } /* end for */

  dstartx = pcam->dxc[pcam->inumcentralpixels - 1] + dhsize;
  dstarty = pcam->dyc[pcam->inumcentralpixels - 1] + dtsize;

  if(pcam->inumgappixels > 0){   /* generate the positions of the gap pixels */
    
    j = pcam->inumcentralpixels;

    for(i=0; i<pcam->inumgappixels; i=i+6){
      pcam->dxc[j + i ] = dstartx + 2. * (i/6 + 1) * dpsize; 
      pcam->dyc[j + i ] = dstarty;
      pcam->dpixsizefactor[j + i] = 1.;
      pcam->dxc[j + i + 1] = pcam->dxc[j + i ] / 2.;
      pcam->dyc[j + i + 1] = sqrt(3.) * pcam->dxc[j + i + 1];
      pcam->dpixsizefactor[j + i + 1] = 1.;
      pcam->dxc[j + i + 2] = - pcam->dxc[j + i + 1];
      pcam->dyc[j + i + 2] = pcam->dyc[j + i + 1];
      pcam->dpixsizefactor[j + i+ 2] = 1.;
      pcam->dxc[j + i + 3] = - pcam->dxc[j + i];
      pcam->dyc[j + i + 3] = dstarty;
      pcam->dpixsizefactor[j + i+ 3] = 1.;
      pcam->dxc[j + i + 4] = pcam->dxc[j + i + 2];
      pcam->dyc[j + i + 4] = - pcam->dyc[j + i + 2];
      pcam->dpixsizefactor[j + i+ 4] = 1.;
      pcam->dxc[j + i + 5] = pcam->dxc[j + i + 1];
      pcam->dyc[j + i + 5] = - pcam->dyc[j + i + 1];
      pcam->dpixsizefactor[j + i + 5] = 1.;
    } /* end for */
  } /* end if */

  /* generate positions of the outer pixels */

  if( pcam->inumbigpixels > 0 ){

    j = pcam->inumcentralpixels + pcam->inumgappixels;

    for(i=0; i<pcam->inumbigpixels; i++){
      pcam->dpixsizefactor[j + i] = 2.;
    }

    in = 0;

    nrow = (int) ceil(dstartx / 2. / dpsize);    

    while(in < pcam->inumbigpixels){

      pcam->dxc[j + in] = dstartx + dpsize;
      pcam->dyc[j + in] = dstarty + 2 * dpsize / sqrt(3.);
      pcam->dxc[j + in + nrow] = dstartx / 2. - dpsize / 2.;
      pcam->dyc[j + in + nrow] = sqrt(3.)/2. * dstartx + 2.5 * dpsize/sqrt(3.);
      pcam->dxc[j + in + 3 * nrow - 1] = - pcam->dxc[j + in];
      pcam->dyc[j + in + 3 * nrow - 1] = pcam->dyc[j + in];
      pcam->dxc[j + in + 3 * nrow] = - pcam->dxc[j + in];
      pcam->dyc[j + in + 3 * nrow] = - pcam->dyc[j + in];
      pcam->dxc[j + in + 5 * nrow - 1] = pcam->dxc[j + in + nrow];
      pcam->dyc[j + in + 5 * nrow - 1] = - pcam->dyc[j + in + nrow];
      pcam->dxc[j + in + 6 * nrow - 1] = pcam->dxc[j + in];
      pcam->dyc[j + in + 6 * nrow - 1] = - pcam->dyc[j + in];
      for(i=1; i<nrow; i++){
	pcam->dxc[j + in + i] = pcam->dxc[j + in] - i * dpsize;
	pcam->dyc[j + in + i] = pcam->dyc[j + in] + i * dpsize * sqrt(3.);
	pcam->dxc[j + in + i + nrow] = pcam->dxc[j + in + nrow] - i * 2 * dpsize;
	pcam->dyc[j + in + i + nrow] = pcam->dyc[j + in + nrow];
        pcam->dxc[j + in + 3 * nrow - 1 - i] = - pcam->dxc[j + in + i];
	pcam->dyc[j + in + 3 * nrow - 1- i] = pcam->dyc[j + in + i];
	pcam->dxc[j + in + i + 3 * nrow] = - pcam->dxc[j + in + i];
	pcam->dyc[j + in + i + 3 * nrow] = - pcam->dyc[j + in + i];
	pcam->dxc[j + in + 5 * nrow - 1 - i] = pcam->dxc[j + in + i + nrow];
	pcam->dyc[j + in + 5 * nrow - 1 - i] = - pcam->dyc[j + in + i + nrow];
	pcam->dxc[j + in + 6 * nrow - 1 - i] = pcam->dxc[j + in + i];
	pcam->dyc[j + in + 6 * nrow - 1 - i] = - pcam->dyc[j + in + i];
      }
      in = in + 6 * nrow;
      dstartx = dstartx + 2. * dpsize;
      nrow = nrow + 1;
    } /* end while */

  } /* 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);

}
//!@}

//!-----------------------------------------------------------
// @name bpoint_is_in_pix
//                                             
// @desc check if a point (x,y) in camera coordinates is inside a given pixel
// 
// @var *pcam     structure camera containing all the
//                camera information
// @var dx, dy    point coordinates in centimeters
// @var ipixnum   pixel number (starting at 0)
// @return        TRUE if the point is inside the pixel, FALSE otherwise
//
// DP
//
// @date Thu Oct 14 16:59:04 CEST 1999
//------------------------------------------------------------
// @function  

//!@{

/******** bpoint_is_in_pix() ***************************************/

int bpoint_is_in_pix(double dx, double dy, int ipixnum, struct camera *pcam){
  /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
  /* the pixel is assumed to be a "closed set" */

  double a, b; /* a = length of one of the edges of one pixel, b = half the width of one pixel */
  double c, xx, yy; /* auxiliary variable */

  b = pcam->dpixdiameter_cm / 2. * pcam->dpixsizefactor[ipixnum];
  a = pcam->dpixdiameter_cm / sqrt(3.) * pcam->dpixsizefactor[ipixnum];
  c = 1. - 1./sqrt(3.);
  if((ipixnum < 0)||(ipixnum >= pcam->inumpixels)){
    fprintf(stderr, "Error in bpoint_is_in_pix: invalid pixel number %d\n", ipixnum);
    fprintf(stderr, "Exiting.\n");
    exit(203);
  }
  xx = dx - pcam->dxc[ipixnum];
  yy = dy - pcam->dyc[ipixnum];

  if(((-b <= xx) && (xx <= 0.) && ((-c * xx - a) <= yy) && (yy <= ( c * xx + a))) ||
     ((0. <  xx) && (xx <= b ) && (( c * xx - a) <= yy) && (yy <= (-c * xx + a)))   ){
    return(TRUE); /* point is inside */
  }
  else{
    return(FALSE); /* point is outside */
  }
}

//!@}

//------------------------------------------------------------
// @name dist_r_P                          
//                                     
// @desc distance straight line r - point P
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//------------------------------------------------------------
// dist_r_P
//
// distance straight line r - point P
//------------------------------------------------------------

float 
dist_r_P(float a, float b, float c, 
         float l, float m, float n,
         float x, float y, float z)
{
  return (
          sqrt((SQR((a-x)*m-(b-y)*l) +
                SQR((b-y)*n-(c-z)*m) +
                SQR((c-z)*l-(a-x)*n))/
               (SQR(l)+SQR(m)+SQR(n))
               )
          );
}
// @endcode


//=------------------------------------------------------------
//!@subsection Log of this file.

//!@{
//
// $Log: not supported by cvs2svn $
// Revision 1.3  1999/10/22 15:01:28  petry
// version sent to H.K. and N.M. on Fri Oct 22 1999
//
// Revision 1.2  1999/10/22 09:44:23  petry
// first synthesized version which compiles and runs without crashing;
//
// Revision 1.1.1.1  1999/10/21 16:35:10  petry
// first synthesised version
//
// Revision 1.13  1999/03/15  14:59:05  gonzalez
// camera-1_1
//
// Revision 1.12  1999/03/02  09:56:10  gonzalez
// *** empty log message ***
//
//
//!@}

//=EOF
