source: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx@ 1223

//!/////////////////////////////////////////////////////////////////////
//
// 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.32 $
// $Author: blanch $ 
// $Date: 2002-02-28 15:04:52 $
//
////////////////////////////////////////////////////////////////////////
// @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 "TRandom.h"
#include "TApplication.h"

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

#include "TArrayC.h"

#include "MTrigger.hxx"
#include "MFadc.hxx"
#include "MLons.hxx"

#include "MRawRunHeader.h"
#include "MRawEvtData.h"
#include "MRawEvtHeader.h"
#include "MRawCrateData.h"
#include "MMcEvt.hxx"
#include "MMcTrig.hxx"
#include "MMcRunHeader.hxx"
#include "MMcTrigHeader.hxx"
#include "MMcFadcHeader.hxx"

/*!@" 

  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__


//!@}

//=-----------------------------------------------------------
//!@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;

static int Write_McEvt  = TRUE;
static int Write_McTrig = TRUE;
static int Write_McFADC = TRUE;
static int Write_RawEvt = 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; 

//@: flag: TRUE: show all fadc singnal in the screen; FALSE: don't
static int FADC_Scan = FALSE;

//@: flag: TRUE: show all trigger singnal in the screen; FALSE: don't
static int Trigger_Scan = FALSE;

//@: flag: TRUE: loop trigger analysis over several thresholds, multiplicities and topologies; FALSE: a single trigger configuration
static int Trigger_Loop = FALSE;

//@: flag: TRUE: Different threshold for each pixel ; FALSE: same threshold for all pixels
static int Individual_Thres_Pixel = FALSE;

//@: Properties of the trigger
static float Trigger_gate_length = 6.0;
static float Trigger_response_ampl = 1.0;
static float Trigger_response_fwhm = 2.0;
static float Trigger_overlaping_time= 0.25;
static float Trigger_noise= 0.3;

//@: Properties of the FADC
static float FADC_response_ampl = MFADC_RESPONSE_AMPLITUDE;
static float FADC_response_fwhm = MFADC_RESPONSE_FWHM;
static float FADC_noise = 2.0;

//@: Trigger conditions for a single trigger mode
static float qThreshold[TRIGGER_PIXELS];  
static int Trigger_multiplicity = 4;
static int Trigger_topology = 2;

//@: Upper and lower edges of the trigger loop
static float Trigger_loop_lthres = 2.0;
static float Trigger_loop_uthres = 10.0;
static float Trigger_loop_sthres = 1.0;
static int Trigger_loop_lmult = 2;
static int Trigger_loop_umult = 10;
static int Trigger_loop_ltop = 0;
static int Trigger_loop_utop = 2;

//!@}

/*!@"

  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


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

 
//!@}

/*!@"

  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|.

  @"*/


//!@}

extern char FileName[];


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

//!@{

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

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

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

  char inname[256];           //@< input file name
  char starfieldname[256];    //@< starfield input file name

  char datname[256];          //@< data (ASCII) output file name

  char rootname[256] ;        //@< ROOT file name 
  char rootname_loop[256] ;   //@< ROOT file name 

  char parname[256];          //@< parameters file name

  char nsbpathname[256];      //@< directory with the dataabse for th NSB

  char flag[SIZE_OF_FLAGS + 1];  //@< flags in the .rfl 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)

  int inumphe;                //@< number of photoelectrons in an event from showers
  float inumphensb;             //@< number of photoelectrons in an event from nsb

  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
  float thetashw, phishw;     //@< parameters of a given shower
  float coreD, coreX, coreY;         //@< core position
  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 photons in a given run
  int ntcph=0;                //@< total number of photons

  int i, j, k;                //@< simple counters

  int addElecNoise;           //@< Will we add ElecNoise?

  int simulateNSB;            //@< Will we simulate NSB?
  int nphe2NSB;               //@< From how many phe will we simulate NSB?
  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 nsb_phepns[iMAXNUMPIX];      //@< NSB values for each pixel

  Float_t nsb_trigresp[TRIGGER_TIME_SLICES];    //@< array to write the trigger
                                                //@< response from the database
  Float_t nsb_fadcresp[(Int_t) SLICES_MFADC];   //@< array to write the fadc
                                                //@< response from the database
  Byte_t trigger_map[((Int_t)(TRIGGER_PIXELS/8))+1];//@< Pixels on when the
                                                    //@< camera triggers
  Float_t fadc_elecnoise[CAMERA_PIXELS];  //@< Electronic niose for each pixel
  Float_t fadc_pedestals[CAMERA_PIXELS];  //@< array for fadc pedestals values

  float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
    EXTWAVEBAND1,
    EXTWAVEBAND2,
    EXTWAVEBAND3,
    EXTWAVEBAND4,
    EXTWAVEBAND5
  };                          
  float zenfactor=1.0;  //  correction factor calculated from the extinction

  float qTailCut;             //@< Tail Cut value
  int nIslandsCut;            //@< Islands Cut value
  int countIslands=FALSE;           //@< Will we count the islands?
  int anaPixels;
    
  float fCorrection;          //@< Factor to apply to pixel values (def. 1.)

  int ntrigger = 0;           //@< number of triggers in the whole file
  int btrigger = 0;           //@< trigger flag
  int ithrescount;            //@< counter for loop over threshold trigger
  float fthrescount;          //@< value for loop over threshold trigger
  int imulticount;            //@< counter for loop over multiplicity trigger
  int itopocount;             //@< counter for loop over topology trigger
  int isorttopo[3];           //@< sorting the topologies
  int icontrigger;            //@< number of trigger conditions to be analised
  UShort_t numPix;            //@< number of sets of fadc written counts

  float fpixelthres[TRIGGER_PIXELS];         //@< Threshold values

  TArrayC *fadcValues;  //@< the analog Fadc signal for pixels

  float plateScale_cm2deg;    //@< plate scale (deg/cm)
  float degTriggerZone;       //@< trigger area in the camera (radius, in deg.)

  int still_in_loop = FALSE;

  //@< Variables to fill the McRunHeader
  Int_t sfRaH = 5;
  Int_t sfRaM = 34;
  Int_t sfRaS = 32;
  Int_t sfDeD = 22;
  Int_t sfDeM = 00;
  Int_t sfDeS = 52;
  Float_t shthetamax = 0.0;
  Float_t shthetamin = 0.0;
  Float_t shphimax = 0.0;
  Float_t shphimin = 0.0;
  Float_t telestheta = 0.0;
  Float_t telesphi = 0.0;
  Float_t sofftheta = 0.0;
  Float_t soffphi = 0.0;
  UInt_t corsika = 5200 ;


  struct camera cam; // structure holding the camera definition


  //!@' @#### Definition of variables for |getopt()|.
  //@'

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

  //!@' @#### Initialisation of some variables
  //@'

  for(i=0;i<iMAXNUMPIX;i++){
    nsb_phepns[i]=0;
    for(j=0;j<iNUMWAVEBANDS;j++)
      nsbrate_phepns[i][j]=0.0;    //@< Starlight rates initialised at 0.0
  }
  /*!@'

    @#### 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();

  Write_McEvt  = get_write_McEvt()  ; 
  Write_McTrig = get_write_McTrig() ; 
  Write_McFADC = get_write_McFadc() ; 
  Write_RawEvt = get_write_RawEvt() ; 

  FADC_Scan = get_FADC_Scan();
  Trigger_Scan = get_Trigger_Scan();

  Individual_Thres_Pixel = get_indi_thres_pixel();

  get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm);

  get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm);

  Trigger_Loop = get_Trigger_Loop(&Trigger_loop_lthres, &Trigger_loop_uthres, &Trigger_loop_sthres, &Trigger_loop_lmult, &Trigger_loop_umult, &Trigger_loop_ltop, &Trigger_loop_utop);

  icontrigger =((int)((Trigger_loop_uthres-Trigger_loop_lthres)
                /Trigger_loop_sthres)+1)*
    (Trigger_loop_umult-Trigger_loop_lmult+1)*
    (Trigger_loop_utop-Trigger_loop_ltop+1);
  
  if (!Trigger_Loop){
    get_Trigger_Single (qThreshold, &Trigger_multiplicity, &Trigger_topology);
    icontrigger=1;
  }
  else
    get_threshold(qThreshold);

  // get filenames

  strcpy( inname, get_input_filename() );
  strcpy( starfieldname, get_starfield_filename() );
  strcpy( datname, get_data_filename() );
  strcpy( rootname, get_root_filename() );
  strcpy( rootname_loop, get_loop_filename() );
  strcpy( ct_filename, get_ct_filename() );
  strcpy( nsbpathname, get_nsb_directory() );

  // get different parameters of the simulation

  qTailCut = get_tail_cut();
  addElecNoise = add_elec_noise(&FADC_noise, &Trigger_noise);
  simulateNSB = get_nsb( &meanNSB, &nphe2NSB );  
  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, 
      "Stars", starfieldname,
      "NSB database", nsbpathname,
      "CT", ct_filename, 
      "Data", datname,
      "ROOT",  rootname 
      );

  // log Trigger information

  if (Trigger_Loop) {
    log(SIGNATURE,
        "%s:\n\t%20s: from %5.2f to %5.2f by %5.2f step\n\t%20s: %i - %i\n\t%20s: %i - %i\n\t%20s\n",
        "Trigger Loop mode",
        "Threshold",Trigger_loop_lthres,Trigger_loop_uthres,Trigger_loop_sthres,
        "Multiplicity",Trigger_loop_lmult,Trigger_loop_umult,
        "Topology",Trigger_loop_ltop,Trigger_loop_utop,
        rootname_loop);
  }
  else if (Individual_Thres_Pixel == FALSE){
    log(SIGNATURE,
        "%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n",
        "Single Trigger mode",
        "Threshold",qThreshold[0],
        "Multiplicity",Trigger_multiplicity,
        "Topology",Trigger_topology);
  }   
  else{
    log(SIGNATURE,
        "%s:\n\t%20s: %s\n\t%20s: %i\n\t%20s: %i\n",
        "Single Trigger mode",
        "Threshold","Different for each pixel",
        "Multiplicity",Trigger_multiplicity,
        "Topology",Trigger_topology);
  }
  // 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 flags information

  log(SIGNATURE,
      "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n",
      "Root output",
      "Write_McEvt",   ONoff(Write_McEvt),  
      "Write_McTrig",  ONoff(Write_McTrig),  
      "Write_McFADC",  ONoff(Write_McFADC),  
      "Write_RawEvt",  ONoff(Write_RawEvt));
      
  // log parameters information
  
  log(SIGNATURE,
      "%s:\n\t%20s: %f\n\t%20s: %f %s\n",
      "Parameters",
      "t0 (Tail-cut)", qTailCut,
      "NSB (phes/pixel)", meanNSB, ONoff(simulateNSB));
  
  // log selections
  
  log(SIGNATURE,
      "%s:\n\t%20s: %s  (%f:%f)\n",
      "Selections:",
      "Energy", ONoff(Select_Energy), Select_Energy_le, Select_Energy_ue);
  
  //  Definition and initialization of array to save trigger statistics
 
  int ***ntriggerloop;

  ntriggerloop= new int ** [(int)((Trigger_loop_uthres-Trigger_loop_lthres)
                           /Trigger_loop_sthres)]; 
  for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;fthrescount+=Trigger_loop_sthres, ithrescount++){
    ntriggerloop[ithrescount]= new int * [Trigger_loop_umult-Trigger_loop_lmult+1];
    for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
      ntriggerloop[ithrescount][imulticount]= new int [Trigger_loop_utop-Trigger_loop_ltop+1];
      for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
        ntriggerloop[ithrescount][imulticount][itopocount]=0;
      }
    }   
  }

  // We should be careful that topologies are sort from 
  // the less to the more restrictive one. 

  if (Trigger_loop_utop==Trigger_loop_ltop)
    for(int is=0; is<3;is++)
      isorttopo[is]=is;
  else {
      isorttopo[0]=1;
      isorttopo[1]=0;
      isorttopo[2]=2;
  }    

  // 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 camera setup

  read_pixels(&cam);

  Int_t Lev0, Lev1, Lev2 ; 
  
  fadcValues =  new TArrayC(FADC_SLICES);

  // number of pixels for parameters
    
  anaPixels = get_ana_pixels();
  anaPixels = (anaPixels == -1) ? ct_NPixels : anaPixels;

  // initialise ROOT

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

  // initialise instance of Trigger and FADC classes

  MTrigger  Trigger(Trigger_gate_length, Trigger_overlaping_time,
                    Trigger_response_ampl, Trigger_response_fwhm);  //@< A instance of the Class MTrigger 

  //  Initialise McTrig information class if we want to save trigger informtion

  MMcTrig **McTrig = NULL; 
  MMcTrigHeader **HeaderTrig = NULL; 
  MMcFadcHeader **HeaderFadc = NULL;


  if (Write_McTrig){

    McTrig = new MMcTrig * [icontrigger];
  
    for (i=0;i<icontrigger;i++) {
      McTrig[i] = new MMcTrig();
    }

    HeaderTrig = new MMcTrigHeader * [icontrigger];
  
    for (i=0;i<icontrigger;i++) {
      HeaderTrig[i] = new MMcTrigHeader();
    }
  }

  if (Write_McFADC){

    HeaderFadc = new MMcFadcHeader * [icontrigger];
  
    for (i=0;i<icontrigger;i++) {
      HeaderFadc[i] = new MMcFadcHeader();
    }
  }

  MFadc fadc(FADC_response_ampl,FADC_response_fwhm) ; //@< A instance of the Class MFadc

  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //
  // Set the FADC pedestals for that run
  //  Some modifications
  // mut be done to simulate a more realistic distribution of the pedestals.
  // This simualtion is done int the SetPedestals methode inside the 
  // class MFadc
  // Currentlly a input_pedestal array is declared with the pedestals.
  // Thy can also be set randomly following a flat distribution.
  //
  /////////////////////////////////////////////////////////////////////

  Float_t input_pedestals[CAMERA_PIXELS];

  for(i=0;i<CAMERA_PIXELS;i++)
    input_pedestals[i]=10.5;
  
  fadc.SetPedestals(input_pedestals);

  // Prepare the raw data output

       // Header Tree

  MRawRunHeader *RunHeader= new MRawRunHeader();
  MMcRunHeader  *McRunHeader = new MMcRunHeader();

       // Header branch

  MRawEvtHeader **EvtHeader = NULL;

  if (Write_RawEvt) {
    EvtHeader = new MRawEvtHeader * [icontrigger]; 

    for (i=0;i<icontrigger;i++) {
      EvtHeader[i] = new MRawEvtHeader();
    }
  }

       // Data branch

  MRawEvtData **EvtData = NULL;    // Data branch

  if (Write_RawEvt) {
    EvtData = new MRawEvtData * [icontrigger]; 

    for (i=0;i<icontrigger;i++) {
      EvtData[i] = new MRawEvtData();
      EvtData[i]->Init(RunHeader);     //  We need thr RunHeader to read
                                       //  number of pixels
    }
  }

  MMcEvt  *McEvt = new MMcEvt (); 

  //  
  // initalize a temporal ROOT file 
  //
  
  TFile outfile_temp ( rootname , "RECREATE" ); 

  //      create a Tree for the Header Event
  TTree HeaderTree("RunHeaders","Headers of Run");
  
  //      define branches of Header Tree

  char help[4];  

  HeaderTree.Branch("MRawRunHeader","MRawRunHeader",
                    &RunHeader);

  HeaderTree.Branch("MMcRunHeader","MMcRunHeader",
                    &McRunHeader);

  if(!Trigger_Loop && Write_McTrig){
    
    HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader", 
                      &HeaderTrig[0]);
  }
  if (Trigger_Loop && Write_McTrig){
    for(char branchname[20],i=0;i<icontrigger;i++){
      
      sprintf(help,"%i",i+1);
      strcpy (branchname, "MMcTrigHeader;");
      strcat (branchname, & help[0]);
      strcat (branchname, ".");
      HeaderTree.Branch(branchname,"MMcTrigHeader", 
                        &HeaderTrig[i]);
    }
  }  

  if(!Trigger_Loop && Write_McFADC){
    
    HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader", 
                      &HeaderFadc[0]);
  }
  if (Trigger_Loop && Write_McFADC){
    for(char branchname[20],i=0;i<icontrigger;i++){
      
      sprintf(help,"%i",i+1);
      strcpy (branchname, "MMcFadcHeader;");
      strcat (branchname, & help[0]);
      strcat (branchname, ".");
      HeaderTree.Branch(branchname,"MMcFadcHeader", 
                        &HeaderFadc[i]);
    }
  }  

  //  Fill branches for MRawRunHeader

  RunHeader->SetMagicNumber(kMagicNumber);
  RunHeader->SetFormatVersion(2);
  RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
  RunHeader->SetRunType(256);
  RunHeader->SetRunNumber(0);
  RunHeader->SetNumSamples(0,FADC_SLICES);
  RunHeader->SetNumCrates(1);
  RunHeader->SetNumPixInCrate(CAMERA_PIXELS);


  //  Fill branches for MMcTrigHeader
  
  if(!Trigger_Loop && Write_McTrig){

    HeaderTrig[0]->SetTopology((Short_t) Trigger_topology);
    HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity);
    HeaderTrig[0]->SetThreshold(qThreshold);
    HeaderTrig[0]->SetAmplitud(Trigger_response_ampl);
    HeaderTrig[0]->SetFwhm(Trigger_response_fwhm);
    HeaderTrig[0]->SetOverlap(Trigger_overlaping_time);
    HeaderTrig[0]->SetGate(Trigger_gate_length);
  }
  if(Trigger_Loop && Write_McTrig){

    int iconcount;
    for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++,fthrescount+=Trigger_loop_sthres){
      for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
        for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
          HeaderTrig[iconcount]->SetTopology((Short_t) isorttopo[itopocount+Trigger_loop_ltop]);
          HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult);
          for(int i=0;i<TRIGGER_PIXELS;i++){
            fpixelthres[i]=
              ((Float_t)(fthrescount)>=qThreshold[i])?
              (Float_t)(fthrescount):qThreshold[i];
              }
          HeaderTrig[iconcount]->SetThreshold( fpixelthres);
          HeaderTrig[iconcount]->SetAmplitud(Trigger_response_ampl);
          HeaderTrig[iconcount]->SetFwhm(Trigger_response_fwhm);
          HeaderTrig[iconcount]->SetOverlap(Trigger_overlaping_time);
          HeaderTrig[iconcount]->SetGate(Trigger_gate_length);
          iconcount++;
        }
      }   
    }
  }

  //  Fill branches for MMcFadcHeader

  for(i=0;i<CAMERA_PIXELS;i++){
      fadc_elecnoise[i]=FADC_noise;
  }

  fadc.GetPedestals(&fadc_pedestals[0]);

  if(!Trigger_Loop && Write_McFADC){

    HeaderFadc[0]->SetShape(0.0);
    HeaderFadc[0]->SetAmplitud(FADC_response_ampl);
    HeaderFadc[0]->SetFwhm(FADC_response_fwhm);
    HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],anaPixels);
    HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0],anaPixels);
  }
  if(Trigger_Loop && Write_McFADC){
    int iconcount;
    for (iconcount=0,ithrescount=0,fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
      for (imulticount=0;imulticount<=Trigger_loop_umult-Trigger_loop_lmult;imulticount++){
        for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
          HeaderFadc[iconcount]->SetShape(0.0);
          HeaderFadc[iconcount]->SetAmplitud(Trigger_response_ampl);
          HeaderFadc[iconcount]->SetFwhm(Trigger_response_fwhm);
          HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0],anaPixels);
          HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],anaPixels);
          iconcount++;
        }
      }   
    }
  }

  //  Fill the Header Tree with the current leaves of each branch
  // HeaderTree.Fill() ;
            

  //      create a Tree for the Event data stream 
  TTree EvtTree("Events","Normal Triggered Events");

  if (Write_McEvt){

    EvtTree.Branch("MMcEvt","MMcEvt", 
                   &McEvt);
  }

  if(!Trigger_Loop){
    
    if (Write_RawEvt){
      EvtTree.Branch("MRawEvtHeader","MRawEvtHeader", 
                     &EvtHeader[0]);
      EvtTree.Branch("MRawEvtData","MRawEvtData", 
                     &EvtData[0]);
    }
    if (Write_McTrig){
      EvtTree.Branch("MMcTrig","MMcTrig", 
                     &McTrig[0]);
    }    
  }
  else{
    if (Write_McTrig){
      for(char branchname[10],i=0;i<icontrigger;i++){
      
        sprintf(help,"%i",i+1);
        strcpy (branchname, "MMcTrig;");
        strcat (branchname, & help[0]);
        strcat (branchname, ".");
        EvtTree.Branch(branchname,"MMcTrig", 
                       &McTrig[i]);
      }
    }
  }  

  if (Trigger_Loop && Write_RawEvt){
    for(char branchname[15],i=0;i<icontrigger;i++){
      
      sprintf(help,"%i",i+1);
      strcpy (branchname, "MRawEvtHeader;");
      strcat (branchname, & help[0]);
      strcat (branchname, ".");
      EvtTree.Branch(branchname,"MRawEvtHeader", 
                     &EvtHeader[i]);
    }
    for(char branchname[15],i=0;i<icontrigger;i++){
      
      sprintf(help,"%i",i+1);
      strcpy (branchname, "MRawEvtData;");
      strcat (branchname, & help[0]);
      strcat (branchname, ".");
      EvtTree.Branch(branchname,"MRawEvtData", 
                     &EvtData[i]);
    }
  }

  TApplication theAppTrigger("App", &argc, argv);

  if(FADC_Scan){
      if (gROOT->IsBatch()) {
          fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
          //    return 1;
      }
  }
  if(FADC_Scan){
      //TApplication theAppFadc("App", &argc, argv);

      if (gROOT->IsBatch()) {
          fprintf(stderr, "%s: cannot run in batch mode\n", argv[0]);
          //    return 1;
      }
  }
  
  // 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;

  // prepare the NSB simulation

  //  Instance of the Mlons class
  MLons lons(Trigger_response_ampl, Trigger_response_fwhm,
             MFADC_RESPONSE_AMPLITUDE, MFADC_RESPONSE_FWHM);

  lons.SetPath(nsbpathname);

  if( simulateNSB ){

    //
    // Calculate the non-diffuse NSB photoelectron rates
    //
    
    k = produce_nsbrates( starfieldname,
                          &cam,
                          nsbrate_phepns);

    if (k != 0){
      cout << "Error when reading starfield... \nExiting.\n";
      exit(1);
    }
  
    // 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];
    }

    // calculate nsb rate including diffuse and starlight
    // we also include the extinction effect
    for(j=0;j<iNUMWAVEBANDS;j++){
        // calculate the effect of the atmospheric extinction 
    
        zenfactor = pow(10., -0.4 * ext[j] );
        
        for(i=0; i<cam.inumpixels;i++){
            nsb_phepns[i]+=diffnsb_phepns[i]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[i][j];
        }
    }
  }

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

  if(check_reflector_file( inputfile )==FALSE){
    exit(1);
  }

  // open data file
  
  log( SIGNATURE, "Opening data \"dat\" file %s\n", datname );
  datafile.open( datname );
  
  if ( datafile.bad() ) 
    error( SIGNATURE, "Cannot open data file: %s\n", datname );
  
  // initializes flag
  
  strcpy( flag, "                                        \0" );

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

     
  //!@' @#### Main loop.
  //@'

  // get flag
    
  fread( flag, SIZE_OF_FLAGS, 1, inputfile );

  // loop over the file

  still_in_loop = TRUE;

  while (
         ((! Data_From_STDIN) && ( !feof(inputfile) ))
         ||
         (Data_From_STDIN && still_in_loop)
         ) { 

    // reading .rfl files 
    if(!isA( flag, FLAG_START_OF_RUN )){

      //  We exit
      //error( SIGNATURE, "Expected start of run flag, but found: %s\n", flag );
      //  We break the main loop
      cout<<"Warning: Expected start of run flag, but found:"<<flag<<endl;
      cout<<"         We break the main loop"<<endl;
      break;      
    }
    else { // found start of run

      nshow=0;

      fread( flag, SIZE_OF_FLAGS, 1, inputfile );

      while( isA( flag, FLAG_START_OF_EVENT   )){ // while there is a next event
        //
        // Clear Trigger and Fadc
        //
        Trigger.Reset() ; 
        Trigger.ClearFirst();
        Trigger.ClearZero();
        fadc.Reset() ; 

        ++nshow;
        log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
        
        // get MCEventHeader
        fread( (char*)&mcevth, mcevth.mysize(), 1, 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;
          
        }

        // energy cut
        
        if ( Select_Energy ) {
          if (( mcevth.get_energy() < Select_Energy_le ) ||
              ( mcevth.get_energy() > Select_Energy_ue )) {
            log(SIGNATURE, "select_energy: shower rejected.\n");
            continue;
          }
        }

        // Read first and last time and put inumphe to 0

        mcevth.get_times(&arrtmin_ns,&arrtmax_ns);
        inumphe=0;

        // read the photons and produce the photoelectrons

        k = produce_phes( inputfile,
                          &cam,
                          WAVEBANDBOUND1,
                          WAVEBANDBOUND6,
                          &Trigger,   // will be changed by the function!
                          &fadc,      // 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);
        }

        // NSB simulation

        if(simulateNSB && nphe2NSB<=inumphe){

          //  Fill trigger and fadc response in the trigger class from the database
          for(i=0;i<cam.inumpixels;i++)
            if(nsb_phepns[i]>0.0){
              k=lons.GetResponse(nsb_phepns[i],0.1,
                                 & nsb_trigresp[0],& nsb_fadcresp[0]);
              if(k==0){
                cout << "Exiting.\n";
                exit(1);
              }
              Trigger.AddNSB(i,nsb_trigresp);
              fadc.AddSignal(i,nsb_fadcresp);
            }
        }// end if(simulateNSB && nphe2NSB<=inumphe) ...
  
        log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
            ncph, ntcph);

        ntcph += ncph;

        // 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;
        }
        
        cout << "Total number of phes: " << inumphe<<" + ";
        inumphensb=0;
        for (i=0;i<cam.inumpixels;i++){
          inumphensb+=nsb_phepns[i]*TOTAL_TRIGGER_TIME;
        }
        cout<<inumphensb<<endl;

        //++++++++++++++++++++++++++++++++++++++++++++++++++
        // 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).
        //--------------------------------------------------
          
        // TRIGGER HERE

        //
        //   now the noise of the electronic 
        //   (preamps, optical transmission,..)  is introduced. 
        //   This is done inside the class MTrigger by the method ElecNoise. 
        //   
        if (addElecNoise){      
            Trigger.ElecNoise(Trigger_noise) ;
            
            fadc.ElecNoise(FADC_noise) ;
        }

        //  We should simulate the AC coupling behaviour:
        //  For the FADC it is only done for the NSB while producinr
        //  the StarResponse database.
        //  For the trigger is done in the Trigger.Diskriminate(), which
        //  is called later (it should be seperated to speed up the program

        //  now a shift in the fadc signal due to the pedestlas is 
        //  intorduced
        //  This is done inside the class MFadc by the method Pedestals
        fadc.Pedestals();

        //   We study several trigger conditons
        if(Trigger_Loop){

          // Set to zero the flag to know if some conditon has triggered
          btrigger=0;
          //  Loop over trigger threshold
          int iconcount;
          for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
            for (i=0;i<TRIGGER_PIXELS;i++)
              fpixelthres[i]=
                ((Float_t)(fthrescount)>=qThreshold[i])?
                (Float_t)(fthrescount):qThreshold[i];
            Trigger.SetThreshold(fpixelthres);

            Trigger.Diskriminate();
            //
            //   look if in all the signals in the trigger signal branch
            //   is a possible Trigger. Therefore we habe to diskriminate all
            //   the simulated analog signals (Method Diskriminate in class
            //   MTrigger). We look simultanously for the moments at which
            //   there are more than TRIGGER_MULTI pixels above the 
            //   CHANNEL_THRESHOLD. 
            //
            
            //  Set trigger flags to zero
            Lev1=Lev2=0;

            //  loop over multiplicity of trigger configuration
            for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
              Trigger.SetMultiplicity(imulticount);
              Trigger.ClearZero();

              Lev0=(Short_t) Trigger.ZeroLevel();
              if (Lev0>0 || Write_All_Images || btrigger){
                
                //  loop over topologies
                for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
                  Lev1=Lev2=0;
                
                  if(itopocount==0 && imulticount>7) continue;
                  //COBB if(itopocount==2 && imulticount<3) continue;
                  // It only makes to look for a different topology
                  // if there are 3 or more N pixels. 
                  if(imulticount<3)
                    Trigger.SetTopology(1);
                  else
                    {
                      // We should be careful that topologies are sort from 
                      // the less to the more restrictive one. 
                      Trigger.SetTopology(isorttopo[itopocount]);
                    }
                  Trigger.ClearFirst();
                  
                  //
                  //   Start the First Level Trigger simulation
                  //
                  if(Lev0!=0)
                      Lev1=Trigger.FirstLevel();
                  if(Lev1>0) {
                    btrigger= 1;
                    ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]++;
                  }
                  
                  Lev0=1;
                  if(Lev1==0 && (Write_All_Images || btrigger)){
                    btrigger= 1;
                    Lev1=1;
                    Lev0=0;
                  }
                  numPix=0;
                  for (Int_t ii=0;ii<Lev1;ii++){
                      if (Write_McTrig){
                          McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
                          McTrig[iconcount]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
                          Trigger.GetMapDiskriminator(trigger_map);
                          McTrig[iconcount]->SetMapPixels(trigger_map,ii);
                      }
                    //
                    //  fill inside the class fadc the member output
                    //

                    fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
            
                    if( Write_RawEvt ){
                      //
                      //  Fill the header of this event 
                      //
                      
                      EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0);
                      //   fill pixel information
                      for(int i=0;i<ct_NPixels;i++){
                        for (j=0;j<FADC_SLICES;j++){
                          fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
                        }
                        EvtData[iconcount]->AddPixel(i,fadcValues,0);
                      }
                    }
                  }
                  //
                  // Increase counter of analised trigger conditions
                  //
                  iconcount++;
                }
              }
              else{
                break;
              }
            }
            if (!btrigger) break;
          }
          if (btrigger){

            //
            //   fill the MMcEvt with all information  
            //

            if (Write_McEvt) {
              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(),
                           impactD,
                           (UInt_t)mcevth.get_CORSIKA(), 
                           (UInt_t)mcevth.get_AtmAbs(), 
                           (UInt_t)mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot(), 
                           (UInt_t) ncph, 
                           (UInt_t) inumphe,
                           (UInt_t) inumphensb+inumphe);
                }
            //  Fill the Tree with the current leaves of each branch
            i=EvtTree.Fill() ;

            //  Clear the branches
            if(Write_McTrig){
              for(i=0;i<icontrigger;i++){
                McTrig[i]->Clear() ;
              }
            }
            if( Write_RawEvt ){
              for(i=0;i<icontrigger;i++){
                EvtHeader[i]->Clear() ;
                EvtData[i]->DeletePixels();
              }
            }
            if (Write_McEvt)
              McEvt->Clear() ;  
          }
        }
        //  We study a single trigger condition
        else {

          //  Setting trigger conditions
          Trigger.SetMultiplicity(Trigger_multiplicity);
          Trigger.SetTopology(Trigger_topology);
          for (int i=0;i<TRIGGER_PIXELS;i++)
            fpixelthres[i]=qThreshold[i];
          Trigger.SetThreshold(fpixelthres);
                                 
          Trigger.Diskriminate() ; 

          //
          //   look if in all the signals in the trigger signal branch
          //   is a possible Trigger. Therefore we habe to diskriminate all
          //   the simulated analog signals (Method Diskriminate in class
          //   MTrigger). We look simultanously for the moments at which
          //   there are more than TRIGGER_MULTI pixels above the 
          //   CHANNEL_THRESHOLD. 
          //
          
          Lev0 = (Short_t) Trigger.ZeroLevel() ; 
          
          Lev1 = Lev2 = 0 ; 
          
          //
          //   Start the First Level Trigger simulation
          //
          
          if ( Lev0 > 0 || Write_All_Images) {
            Lev1= Trigger.FirstLevel();
          }
          if (Lev1>0){
            ++ntrigger;
          }
          Lev0=1;
          if (Lev1==0 && Write_All_Images){ 
            Lev1=1;
            Lev0=0;
          }

          numPix=0;
          for(Int_t ii=0;ii<Lev1;ii++){
            //  Loop over different level one triggers

            //
            //   fill inside class fadc the member output
            //
            fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));

            if (Write_McTrig){
              McTrig[0]->SetFirstLevel ((ii+1)*Lev0);
              McTrig[0]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
              Trigger.GetMapDiskriminator(trigger_map);
              McTrig[0]->SetMapPixels(trigger_map,ii);
            }

            //  Fill Evt information

            if (Write_RawEvt){

              //
              //  Fill the header of this event 
              //
              
              EvtHeader[0]->FillHeader ( (UShort_t) (ntshow + nshow) ,  20 ) ; 
              
              //   fill pixel information
              
              for(int i=0;i<ct_NPixels;i++){
                for (j=0;j<FADC_SLICES;j++){
                  fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
                }
                EvtData[0]->AddPixel(i,fadcValues,0);
              }
            }       
            //
            //   fill the MMcEvt with all information  
            //
            
            if (Write_McEvt){
              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(),
                           impactD,
                           (UInt_t)mcevth.get_CORSIKA(), 
                           (UInt_t)mcevth.get_AtmAbs(), 
                           (UInt_t)mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot(), 
                           (UInt_t) ncph, 
                           (UInt_t) inumphe,
                           (UInt_t) inumphensb+inumphe);
                }
            //   We don not count photons out of the camera.    
            
            
            //
            //    write it out to the file outfile
            // 
            
            EvtTree.Fill() ; 
                          
            //
            //    if a first level trigger occurred, then 
            //       1. do some other stuff (not implemented) 
            //       2. start the gui tool
            
            if(FADC_Scan){
              if ( Lev0 > 0 ) {
                fadc.ShowSignal( McEvt, (Float_t) 60. ) ; 
              }
            }
            
            if(Trigger_Scan){
              if ( Lev0 > 0 ) {
                Trigger.ShowSignal(McEvt) ;
              }
            }
            
            //    clear all
            if (Write_RawEvt) EvtHeader[0]->Clear() ;
            if (Write_RawEvt) EvtData[0]->DeletePixels();
            if (Write_McEvt) McEvt->Clear() ; 
          } 
          if (Write_McTrig) McTrig[0]->Clear() ;
        }
                
#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 (int 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__
                          
        
        // look for the next event
        
        fread( flag, SIZE_OF_FLAGS, 1, inputfile );

      } // 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);
        
        fread( flag, SIZE_OF_FLAGS, 1, inputfile );
        
        if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
          log(SIGNATURE, "End of file . . .\n");
          still_in_loop  = FALSE;
  
          if ((! Data_From_STDIN) && ( !feof(inputfile) )){
            
            // we have concatenated input files.
            // get signature of the next part and check it.

            if(check_reflector_file( inputfile )==FALSE){
              exit(1);
            }
            
          }
          
          fread( flag, SIZE_OF_FLAGS, 1, inputfile );

        } // end if found end of file
      } // end if found end of run

    } // end if else found start of run
  } // end big while loop

  //<@ Finally we should fill th McRunHeader


  get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
  get_teles_axis(&telestheta,&telesphi);
  get_source_off(&sofftheta,&soffphi);
  mcevth.get_theta_range(&shthetamin, &shthetamax);
  mcevth.get_phi_range(&shphimin,&shphimax);
  corsika=get_corsika_ver();
  if(!Trigger_Loop) icontrigger=0;
  McRunHeader->Fill(icontrigger,
                   !Write_All_Images,
                   Write_McEvt,
                   Write_McTrig,
                   Write_McFADC,
                   CAMERA_PIXELS,
                   (UInt_t)ntshow,
                   (UInt_t)ntrigger,
                   sfRaH,
                   sfRaM,
                   sfRaS,
                   sfDeD,
                   sfDeM,
                   sfDeS,
                   meanNSB,
                   telestheta,
                   telesphi,
                   sofftheta,
                   soffphi,
                   shthetamax,
                   shthetamin,
                   shphimax,
                   shphimin,
                   corsika,
                   (UInt_t)(REFL_VERSION*100),
                   (UInt_t)(VERSION*100));

  //  Fill the Header Tree with the current leaves of each branch
  HeaderTree.Fill() ;
  
  //++
  // put the Event to the root file
  //--

  outfile_temp.Write() ;
  outfile_temp.Close() ; 
  

  // close input file
  
  log( SIGNATURE, "%d event(s), with a total of %d C.photons\n", 
       ntshow, ntcph );
  datafile<<ntshow<<" event(s), with a total of "<<ntcph<<" C.photons"<<endl;
  if (Trigger_Loop){
    log( SIGNATURE, "Trigger Mode. \n");
    log( SIGNATURE, "Fraction of triggers: \n");
    datafile<<"Fraction of triggers: "<<endl;
    for (ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
      for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
        for(itopocount=Trigger_loop_ltop;itopocount<=Trigger_loop_utop;itopocount++){
          log( SIGNATURE, "Thres %5.1f, Multi %d, Topo %d: %5.1f%% (%d out of %d)\n", 
               fthrescount,imulticount,isorttopo[itopocount],((float)ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop] / ((float)ntshow) * 100.0), ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop], ntshow);
          datafile<<"Thres "<<fthrescount<<", Multi "<<imulticount<<", Topo"<<isorttopo[itopocount]<<": ";
          datafile<<((float)ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop] / ((float)ntshow) * 100.0)<<"% ("<<ntriggerloop[ithrescount][imulticount-Trigger_loop_lmult][itopocount-Trigger_loop_ltop]<<" out of "<<ntshow<<")"<<endl;
        }
      }   
    }
  }
  else{
    log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n", 
         ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
    datafile<<"Fraction of triggers: "<<((float)ntrigger) / ((float)ntshow) * 100.0<<" ("<<ntrigger<<" out of "<<ntshow<<" )"<<endl;
  }

  // close files
  
  log( SIGNATURE, "Closing files\n" );
  
  if( ! Data_From_STDIN ){
    fclose( inputfile );
  }
  datafile.close();

  // 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(stdout, "ERROR in %s: ", funct);

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

  perror(funct);

  exit(1);
}
//!@}


//!-----------------------------------------------------------
// @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, icount;
  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

  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);


  // 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", "Opening 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;
  icount = 0;

  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)
    
    if( sscanf(line, "%d %d %f", &i, &j, &qe) != 3 ) 
      break;

    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;
    }

    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

  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] = 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 */

  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./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
  }
}
*/

#define sqrt13 0.577350269 // = 1./sqrt(3.)
#define sqrt3  1.732050807 // = sqrt(3.)

int bpoint_is_in_pix(double dx, double dy, 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" */

    /*
     a = length of one of the edges of one pixel,
     b = half the width of one pixel
     */

    const int    num  = pcam->inumpixels;
    const double diam = pcam->dpixdiameter_cm;

    const double diam2 = diam/2;
    const double diam3 = diam/sqrt3;

    for (int i=0; i<num; i++)
    {
        const double size = pcam->dpixsizefactor[i];

        const double b = diam2 * size;
        const double a = diam3 * size;

        const double xx = dx - pcam->dxc[i];
        const double yy = dy - pcam->dyc[i];

        if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
           ((0. <  xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a)))   ){

            return i; // inside i
        }

        //   return -1; // outside
    }
    return -1; // 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))
               )
          );
}

//------------------------------------------------------------
// @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 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 MTrigger *trigger, // the generated phes
                  class MFadc *fadc,
                  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;
    
  }

  *incph = 0;

  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
  
  
  // loop over the photons
  
  while (1) {

    fread ( flag, SIZE_OF_FLAGS, 1, sp );

    if (isA( flag, FLAG_END_OF_EVENT ))
        break;
          
    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)++;

    //  Chceck if photon is inside trigger time range

    t = photon.get_t() ; 

    if (t-*tmin_ns>TOTAL_TRIGGER_TIME) {
        continue;
    }
          
    //
    // 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 ) ){ 
      continue;
            
    }

    ipixnum = bpoint_is_in_pix(cx, cy, cam);

    // -1 = the photon is in none of the pixels
    //  0 = the phton is in the central pixel, which is not used for trigger
    if (ipixnum==-1 || ipixnum==0) {
        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])){
      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 ) {
      continue;
    }
          
    //+++
    // The photon has produced a photo electron
    //---

    // cout << " accepted\n";
          
    // increment the number of photoelectrons in the relevant pixel
          
    nphe[ipixnum] += 1.0;
          
    // store the new photoelectron

    fadc->Fill(ipixnum,(t-*tmin_ns) , trigger->FillShow(ipixnum,t-*tmin_ns));
    
    *itotnphe += 1;
  }
          
  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
                      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

  MTrigger trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
  MFadc flashadc;
                      
  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" );

  // check if the satrfield input file exists

  if ( infile == NULL ) {

    log( SIGNATURE, "Cannot open starfield input file: %s\n", iname );

    log( SIGNATURE, "There is not NSB from the Stars\n");

    return (0);
  }
  
  // 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],
                              &trigger,  // this is a dummy here
                              &flashadc, // 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],
                      MTrigger *trigger,
                      MFadc *fadc,
                      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

  float t;
  TRandom random;             // Random numbers generator

  random.SetSeed ((UInt_t) (RandomNumber*1000));

  ii = *inphe; // avoid dereferencing
                
  // check if the arrival times are set; if not generate them

  if(*atmin_ns <SIMTIMEOFFSET_NS || *atmin_ns > *atmax_ns){
    *atmin_ns = 0.;
    *atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;

  }
  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 and not too long

    if(simtime_ns< TOTAL_TRIGGER_TIME){
      *atmin_ns = *atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
      *atmax_ns = *atmin_ns + TOTAL_TRIGGER_TIME;
      simtime_ns = TOTAL_TRIGGER_TIME;
    }

    if(simtime_ns> TOTAL_TRIGGER_TIME){
      *atmax_ns =*atmin_ns + TOTAL_TRIGGER_TIME;
      simtime_ns = TOTAL_TRIGGER_TIME;
    }
 
  }

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

      if (inumnsbphe>0.0){
        inumnsbphe = random.Poisson (inumnsbphe );
      }

      // create the photoelectrons
      
      for(k=0; k < inumnsbphe; k++){
        
        t=(RandomNumber * simtime_ns);

        (*fadc).Fill(j,t ,(*trigger).FillNSB(j,t));
        
        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


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

//!@{
//
// $Log: not supported by cvs2svn $
// Revision 1.31  2002/01/18 17:41:02  blanch
// The option of adding noise to all pixels or to not adding the noise
// has been added.
// We removed the pixels larger than 577. When there were more than one
// trigger in one shower, the pixel number was increasing. Now it is
// flagged by the variable MMcTrig::fFirstLvlTrig.
//
// Revision 1.30  2001/11/27 09:49:54  blanch
// Fixing bug which was treating wrongly the extension of star photons.
//
// Revision 1.29  2001/11/14 17:38:23  blanch
// Sveral changes have been done:
//      - bpoint_in_in_pixel has been dodified to speed up the program
//      - Some minor changes have been done to remove warnings
//      - buffer size and split version of the Branches have been removed
//      - Some modifications were needed to adat the program to the new
//        MRawEvtData::DeletePixels
//
// Revision 1.28  2001/10/26 16:31:45  blanch
// Removing several warnings.
//
// Revision 1.27  2001/09/05 10:04:33  blanch
// *** empty log message ***
//
// Revision 1.26  2001/07/19 09:29:53  blanch
// Different threshold for each pixel can be used.
//
// Revision 1.25  2001/05/08 08:07:54  blanch
// New numbering for branches from different trigger conditions has been
// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
// information already filled.
//
// Revision 1.24  2001/04/06 16:48:09  magicsol
// New camera version able to read the new format of the reflector output:
// reflector 0.4
//
// Revision 1.23  2001/03/28 16:13:41  blanch
// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
// been solved.
//
// Revision 1.22  2001/03/20 18:52:43  blanch
// *** empty log message ***
//
// Revision 1.21  2001/03/19 19:53:03  blanch
// Some print outs have been removed.
//
// Revision 1.20  2001/03/19 19:30:06  magicsol
// Minor changes have been done to improve the FADC pedestals treatment.
//
// Revision 1.19  2001/03/05 11:14:41  magicsol
// I changed the position of readinf a parameter. It is a minnor change.
//
// Revision 1.18  2001/03/05 10:36:52  blanch
// A branch with information about the FADC simulation (MMcFadcHeader) is writen
// in the McHeaders tree of the aoutput root file.
// The NSB contribution is only applied if the the number of phe form the shower
// are above a given number.
//
// Revision 1.17  2001/02/23 11:05:57  magicsol
// Small changes due to slightly different output format and the introduction of
// pedesals for teh FADC.
//
// Revision 1.16  2001/01/18 18:44:40  magicsol
// *** empty log message ***
//
// Revision 1.15  2001/01/17 09:32:27  harald
// The changes are neccessary to have the same name for trees in MC and in
// data. So now there should be no differences in MC data and real data from
// FADC system.
//
// Revision 1.14  2001/01/15 12:33:34  magicsol
// Some modifications have been done to use the new (Dec'2000) Raw data format.
// There are also some minnors modifications to adapt the improvements in the
// MTrigger class (overlaping time and trigger cells).
//
// Revision 1.13  2000/10/25 08:14:23  magicsol
// The routine that produce poisson random numbers to decide how many phe
// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
// is used.
//
// Revision 1.12  2000/09/22 17:40:18  harald
// Added a lot of changes done by oscar.
//
// Revision 1.11  2000/09/21 11:47:33  harald
// Oscar found some smaller errors in the calculation of the pixel shape and
// corrected it.
//
// $Log: not supported by cvs2svn $
// Revision 1.31  2002/01/18 17:41:02  blanch
// The option of adding noise to all pixels or to not adding the noise
// has been added.
// We removed the pixels larger than 577. When there were more than one
// trigger in one shower, the pixel number was increasing. Now it is
// flagged by the variable MMcTrig::fFirstLvlTrig.
//
// Revision 1.30  2001/11/27 09:49:54  blanch
// Fixing bug which was treating wrongly the extension of star photons.
//
// Revision 1.29  2001/11/14 17:38:23  blanch
// Sveral changes have been done:
//      - bpoint_in_in_pixel has been dodified to speed up the program
//      - Some minor changes have been done to remove warnings
//      - buffer size and split version of the Branches have been removed
//      - Some modifications were needed to adat the program to the new
//        MRawEvtData::DeletePixels
//
// Revision 1.28  2001/10/26 16:31:45  blanch
// Removing several warnings.
//
// Revision 1.27  2001/09/05 10:04:33  blanch
// *** empty log message ***
//
// Revision 1.26  2001/07/19 09:29:53  blanch
// Different threshold for each pixel can be used.
//
// Revision 1.25  2001/05/08 08:07:54  blanch
// New numbering for branches from different trigger conditions has been
// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
// information already filled.
//
// Revision 1.24  2001/04/06 16:48:09  magicsol
// New camera version able to read the new format of the reflector output:
// reflector 0.4
//
// Revision 1.23  2001/03/28 16:13:41  blanch
// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
// been solved.
//
// Revision 1.22  2001/03/20 18:52:43  blanch
// *** empty log message ***
//
// Revision 1.21  2001/03/19 19:53:03  blanch
// Some print outs have been removed.
//
// Revision 1.20  2001/03/19 19:30:06  magicsol
// Minor changes have been done to improve the FADC pedestals treatment.
//
// Revision 1.19  2001/03/05 11:14:41  magicsol
// I changed the position of readinf a parameter. It is a minnor change.
//
// Revision 1.18  2001/03/05 10:36:52  blanch
// A branch with information about the FADC simulation (MMcFadcHeader) is writen
// in the McHeaders tree of the aoutput root file.
// The NSB contribution is only applied if the the number of phe form the shower
// are above a given number.
//
// Revision 1.17  2001/02/23 11:05:57  magicsol
// Small changes due to slightly different output format and the introduction of
// pedesals for teh FADC.
//
// Revision 1.16  2001/01/18 18:44:40  magicsol
// *** empty log message ***
//
// Revision 1.15  2001/01/17 09:32:27  harald
// The changes are neccessary to have the same name for trees in MC and in
// data. So now there should be no differences in MC data and real data from
// FADC system.
//
// Revision 1.14  2001/01/15 12:33:34  magicsol
// Some modifications have been done to use the new (Dec'2000) Raw data format.
// There are also some minnors modifications to adapt the improvements in the
// MTrigger class (overlaping time and trigger cells).
//
// Revision 1.13  2000/10/25 08:14:23  magicsol
// The routine that produce poisson random numbers to decide how many phe
// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
// is used.
//
// Revision 1.10  2000/07/04 14:10:20  MagicSol
// Some changes have been done in the root output file. The RawEvt tree is only
// stored in the single trigger mode.
// The trigger input parameters are also given by the general input card.
// The diffuse NSB and the star NSB have been decoupled. Now the contribution of
// each one can be studied seperately.
//
// Revision 1.9  2000/06/13 13:25:24  blanch
// The multiple arrays have been replaced, since they do not work
// in alpha machines. Now we are using pointers and the command new
// to allocate memory.
//
// Revision 1.8  2000/05/11 13:57:27  blanch
// The option to loop over several trigger configurations has been included.
// Some non-sense with triggertime range has been solved.
// Montecarlo information and ADC counts are saved in a root file.
// There was a problem with the maximum number of phe that the program could analyse. Now there is not limit.
//
// Revision 1.7  2000/03/24 18:10:46  blanch
// A first FADC simulation and a trigger simulation are already implemented.
// The calculation of the Hillas Parameters have been removed, since it was decided that it should be in the analysis software.
// A loop over trigger threshold and some corretcions in the time range where it looks for a trigger will be implemented as soon as possible.
//
// Revision 1.6  2000/03/20 18:35:11  blanch
// The trigger is already implemented but it does not save the trigger information in any file as it is implemented in timecam. In the next days there will be a version which also creates the files with the trigger information. It is going to be a mixing of the current camera and timecam programs.
//
// Revision 1.5  2000/02/18 17:40:35  petry
// This version includes drastic changes compared to camera.cxx 1.4.
// It is not yet finished and not immediately useful because the
// trigger simulation is not yet re-implemented. I had to take it
// out together with some other stuff in order to tidy the whole
// program up. This is not meant as an insult to anyone. I needed
// to do this in order to be able to work on it.
//
// This version has been put in the repository in order to be
// able to share the further development with others.
//
// If you need something working, wait or take an earlier one.
// See file README.
//
// 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
// of invalid wavelength. This bug is now fixed and the range
// of valid wavelengths extended to 290 - 800 nm.
// This is in preparation for the NSB simulation to come.
// Dirk
//
// Revision 1.2  1999/11/19 08:40:42  harald
// Now it is possible to compile the camera programm under osf1.
//
// Revision 1.1.1.1  1999/11/05 11:59:31  harald
// This the starting point for CVS controlled further developments of the
// camera program. The program was originally written by Jose Carlos. 
// But here you can find a "rootified" version to the program. This means 
// that there is no hbook stuff in it now. Also the output of the
// program changed to the MagicRawDataFormat. 
//
// The "rootification" was done by Dirk Petry and Harald Kornmayer. 
//
// 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