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

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

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

//=-----------------------------------------------------------
//!@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 "TObjArray.h"
                                               
#include "MTrigger.hxx"
#include "MFadc.hxx"
#include "MLons.hxx"

#include "MRawRunHeader.h"
#include "MRawEvtData.h"
#include "MRawEvtHeader.h"
#include "MRawCrateArray.h"
#include "MRawCrateData.h"
#include "MMcEvt.hxx"
#include "MMcEvtBasic.h"
#include "MMcTrig.hxx"
#include "MMcRunHeader.hxx"
#include "MMcConfigRunHeader.h"
#include "MMcCorsikaRunHeader.h"
#include "MMcTrigHeader.hxx"
#include "MMcFadcHeader.hxx"
#include "MGeomCamMagic.h"
#include "MGeomCamMagic919.h"
#include "MGeomCamMagicHG.h"
#include "MGeomCamECO1000.h" 
#include "MGeomCamECO1000HG.h" 
#include "MGeomCamCT1.h"
#include "MGeomCamCT1Daniel.h"
#include "MGeomPix.h"
#include "MGeomCorsikaCT.h"
#include "MTriggerPattern.h"

/*!@" 

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

  @"*/

#include "camera.h"

//!@}

/*!@"

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

  @"*/

//!@{

// flag for debugging (default: OFF )
#define __DEBUG__
#undef  __DEBUG__


//!@}

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

  @"*/

/*!@"

  And this is the information about the whole telescope.

  @"*/

//!@{

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

//@: Number of pixels
static int ct_NPixels;        

//@: Number of CT
static int ct_Number;        

//@: 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 the event headers of all events to output; 
//@:       FALSE: only triggered showers
static int Write_All_Event_Headers = TRUE;

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 signal 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 Int_t FADC_shape = 0;
static float FADC_response_integ = MFADC_RESPONSE_INTEGRAL;
static float FADC_response_fwhm = MFADC_RESPONSE_FWHM;
static Int_t FADC_shape_out = 0;
static float FADC_resp_integ_out = MFADC_RESPONSE_INTEGRAL;
static float FADC_resp_fwhm_out = MFADC_RESPONSE_FWHM;
static float FADC_slices_per_ns = FADC_SLICES_PER_NSEC;
static Int_t FADC_slices_written = FADC_SLICES;
static float FADC_noise_inner = 2.0;
static float FADC_noise_outer = 2.0;
static float DIGITAL_noise = 0.0;
static float FADC_high2low = HIGH2LOWGAIN;

//@: Trigger conditions for a single trigger mode
static float **qThreshold;  
static int Trigger_multiplicity[MAX_NUMBER_OF_CTS];
static int Trigger_topology[MAX_NUMBER_OF_CTS];

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

//@: Direction of each shower
static  float Zenith = 0.0;
static  float Azimutal = 90.0;

//@: Mispointing Simulation
static   int missPointing = 0;
static   float missP_x = 0.0;
static   float missP_y = 0.0;

//@: Point Spread Function Added
static  float Spot_x=0.0;
static  float Spot_y=0.0;
static  float Spotsigma=0.0;

//@: PMT time jitter
static  float pmt_jitter=0.25;

  
//!@}

/*!@"

  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


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

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

//@: table of lightguide = WC;  WC_outer for outer pixels.
static float **WC;
static float **WC_outer;

//@: number of datapoints for the WC curve
static int pointsWC;

//!@}

/*!@"

  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.

  @"*/

/*!@"

  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[];
// switch on  starfield rotation
static int Starfield_rotate = TRUE;

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

//!@{

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

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

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

  char **inname_CT;   //@< array of names for each CT input file
  char starfieldname[256];    //@< starfield input file name
  char qe_filename[256];       //@< qe 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 nsbpath_outer[256];    //@< directory with the dataabse for outer pixels

  char flag[SIZE_OF_FLAGS + 1];  //@< flags in the .rfl file
  char flag_new[SIZE_OF_FLAGS + 1];  //@< New flag for the run header in the .rfl file
  char **GeometryName;   //@< Name of MGeomCam for each CT
  int GeometryCamera[MAX_NUMBER_OF_CTS];  //@< Identification of MGeomCam for each CT
  int TriggerPixels[MAX_NUMBER_OF_CTS]; //@< Number of pixels in the trigger region for each CT

  int reflector_file_version=0;  //@< vrsion of he reflector file

  FILE *inputfile[MAX_NUMBER_OF_CTS]; //@< stream for the input file

  ofstream datafile;          //@< stream for the data file

  MCRunHeader mcrunh;                //@< Run Header class (MC)
  MCEventHeader mcevth[MAX_NUMBER_OF_CTS]; //@< Event Header class (MC)
  MCEventHeader_2 mcevth_2[MAX_NUMBER_OF_CTS];   //@< Event Header class (MC) for reflector > 0.6
  MCCphoton cphoton;          //@< Cherenkov Photon class (MC)

  int inumphe;                //@< number of photoelectrons in an event from showers
  int inumphe_CT[MAX_NUMBER_OF_CTS]; //@< number of photoelectrons in an event from showers
  float inumphensb[MAX_NUMBER_OF_CTS];  //@< 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[MAX_NUMBER_OF_CTS], phiCT[MAX_NUMBER_OF_CTS];  //@< theta and phi of telescopes

  //@: Coordinates of telescopes in Corsika's coordinate system
  float CTx[MAX_NUMBER_OF_CTS];
  float CTy[MAX_NUMBER_OF_CTS];
  float CTz[MAX_NUMBER_OF_CTS];

  float mirror_frac[MAX_NUMBER_OF_CTS];

  float thetashw, phishw;     //@< parameters of a given shower
  float coreX, coreY;         //@< core position
  float impactD[MAX_NUMBER_OF_CTS];  //@< impact parameter
  float l1, m1, n1;           //@< auxiliary variables
  float factorqe_NSB[MAX_NUMBER_OF_CTS]; //@< factor on the NSB depending of QE

  int nshow=0;                //@< partial number of shower in a given run
  int ntshow=0;               //@< total number of showers
  int ncph[MAX_NUMBER_OF_CTS];  //@< partial number of photons in a given run
  int ntcph[MAX_NUMBER_OF_CTS]; //@< total number of photons

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

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

  float riseDiskThres;
  float secureDiskThres;

  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;  
  //@< diffuse NSB values for each pixel derived from meanNSB

  float **nsbrate_phepns;  //@< non-diffuse nsb 
                                                    //@< photoelectron rates 
  float **nsb_phepns;
  //@< NSB values for each pixel

  float **nsb_phepns_rotated;
  //@< NSB values for each pixel after rotation

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

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

  int nstoredevents = 0;
  int flagstoring = 0;

  int ntrigger[MAX_NUMBER_OF_CTS];  //@< 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

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

  TArrayC *fadcValues;     //@< the analog Fadc High gain signal for pixels
  TArrayC *fadcValuesLow;  //@< the analog Fadc Low gain signal for pixels

  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;
  UInt_t corsika = 5200 ;
  Float_t maxpimpact = 0.0;

  // Star Field Rotation variables
  Float_t zenith = 0.0;
  Float_t azimutal = 90.0;
  Float_t rho  , C3 , C2 , C1; 

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

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

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

  inname_CT = new char *[MAX_NUMBER_OF_CTS];
  GeometryName = new char *[MAX_NUMBER_OF_CTS];
  diffnsb_phepns = new float *[MAX_NUMBER_OF_CTS];
  nsb_phepns = new float *[MAX_NUMBER_OF_CTS];
  nsb_phepns_rotated = new float *[MAX_NUMBER_OF_CTS];

  for(int i=0;i<MAX_NUMBER_OF_CTS;i++)
    {
      inname_CT[i] = new char[256];
      GeometryName[i] = new char[50];
      diffnsb_phepns[i] = new float[iMAXNUMPIX];
      nsb_phepns[i] = new float[iMAXNUMPIX]; 
      nsb_phepns_rotated[i] = new float[iMAXNUMPIX]; 
      CTx[i] = 0.; CTy[i] = 0.; CTz[i] = 0.;
    }

  nsbrate_phepns = new float *[iMAXNUMPIX];
  for(int i=0;i<iMAXNUMPIX;i++)
    nsbrate_phepns[i] = new float[iNUMWAVEBANDS];
   
  qThreshold = new float *[MAX_NUMBER_OF_CTS];
  for(int i=0;i<MAX_NUMBER_OF_CTS;i++)
    qThreshold[i] = new float [CAMERA_PIXELS]; 

  for(int i=0;i<iMAXNUMPIX;i++){
    for(int ict=0;ict<MAX_NUMBER_OF_CTS;ict++){      
      nsb_phepns[ict][i]=0;
      nsb_phepns_rotated[ict][i]=0;
    }
    for(j=0;j<iNUMWAVEBANDS;j++)
      nsbrate_phepns[i][j]=0.0;    //@< Starlight rates initialised at 0.0
  }
  for(int i=0;i<MAX_NUMBER_OF_CTS;i++)
    ntcph[i]=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();

  // get number of telescopes and camera geometries

  ct_Number=get_ct_number();
  if (ct_Number > MAX_NUMBER_OF_CTS)
    {
      error( SIGNATURE, "Number of telescopes is larger than maximum allowed (%i). Stoping camera program ...", MAX_NUMBER_OF_CTS);
      exit(1);
    }

  for(int ict=0;ict<ct_Number;ict++){
    ntrigger[ict]=0;
    GeometryCamera[ict]=get_ct_geometry(ict);

    //
    // Get telescope coordinates (if supplied) from  input card (units cm).
    //
    CTx[ict] = get_telescope_location_cm(ict,0);
    CTy[ict] = get_telescope_location_cm(ict,1);
    CTz[ict] = get_telescope_location_cm(ict,2);

    //
    // Get fraction of operative mirror:
    //
    mirror_frac[ict] = get_mirror_fraction(ict);
  }

  // structure holding the camera definition

  TObjArray camgeom;

  for (int i=0; i<ct_Number;i++){
    switch(GeometryCamera[i]){
    case 1: camgeom[i]=new MGeomCamMagic;
      strcpy(GeometryName[i],"MGeomCamMagic");
      TriggerPixels[i]=TRIGGER_PIXELS_1;
      break;
    case 2: camgeom[i]=new MGeomCamMagic919;
      strcpy(GeometryName[i],"MGeomCamMagic919");
      TriggerPixels[i]=TRIGGER_PIXELS_2;
      break;
    case 3: camgeom[i]=new MGeomCamMagicHG;
      strcpy(GeometryName[i],"MGeomCamMagicHG");
      TriggerPixels[i]=TRIGGER_PIXELS_3;
      break;
    case 5: camgeom[i]=new MGeomCamECO1000;
      strcpy(GeometryName[i],"MGeomCamECO1000");
      TriggerPixels[i]=TRIGGER_PIXELS_5;
      break;
    case 6: camgeom[i]=new MGeomCamECO1000HG;
      strcpy(GeometryName[i],"MGeomCamECO1000HG");
      TriggerPixels[i]=TRIGGER_PIXELS_6;
      break;
    case 8: camgeom[i]=new MGeomCamCT1;
      strcpy(GeometryName[i],"MGeomCamCT1");
      TriggerPixels[i]=TRIGGER_PIXELS_8;
      break;
    case 9: camgeom[i]=new MGeomCamCT1Daniel;
      strcpy(GeometryName[i],"MGeomCamCT1Daniel");
      TriggerPixels[i]=TRIGGER_PIXELS_9;
      break;
    default:
      error( SIGNATURE, "Camera geometry %i is not defined. Stoping camera program ...", GeometryCamera[i]);
      exit(1);
      break;
    }
  }
  
  ct_NPixels=0;

  for(int ict=0;ict<ct_Number;ict++)
    ct_NPixels=(((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels()>(UInt_t)ct_NPixels)?
      (Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels():
        ct_NPixels;

  // read WC and QE files

  // FIX ME ! 
  // Currently the PMT N of any camera with same average QE has the same QE.

  QE = new float *** [ct_Number];
  for(int ict=0;ict<ct_Number;ict++){
    strcpy( qe_filename, get_qe_filename(ict));
    read_QE(qe_filename,ict);
  }

  read_WC();

  // write all images, even those without trigger?

  Write_All_Event_Headers = get_write_all_event_headers();

  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_secure_threhold(&riseDiskThres,&secureDiskThres);

  get_FADC_properties
    (&FADC_shape, &FADC_response_integ, &FADC_response_fwhm,
     &FADC_shape_out, &FADC_resp_integ_out, &FADC_resp_fwhm_out,
     &FADC_slices_per_ns, &FADC_slices_written);

  // Allocate memory for the NSB contribution to the FADC signal:
  nsb_fadcresp = new Float_t[(Int_t)(FADC_slices_per_ns*TOTAL_TRIGGER_TIME)];

  FADC_high2low=get_High_to_Low();

  // FIXME --- trigger properties may depend on the Camera geometry

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

  // Trigger loop operation is not implemented for Multi telescopes

  if ( Trigger_Loop && ct_Number > 1 ){
    cout<<"ERROR : camera::main : Trigger loop option is not";
    cout<<" implemented for Multi Telescopes option. "<<icontrigger<<
      " "<<ct_Number<<endl;
    exit(1);
  }
  
  if (!Trigger_Loop){
    get_Trigger_Single (qThreshold, 
                        &Trigger_multiplicity[0], 
                        &Trigger_topology[0]);
    icontrigger=1;
  }
  else
    get_threshold(qThreshold[0]);

  // get filenames

  if (! is_calibration_run())
    for(int ict=0;ict<ct_Number;ict++)
      {
        strcpy( inname_CT[ict], get_input_filename(ict) );
        if (strlen(inname_CT[ict]) == 0)
          {
            printf("\nError: missing input file name for CT id = %d. Exiting.\n\n", ict);
            exit(1);
          }
      }

  strcpy( starfieldname, get_starfield_filename() );
  strcpy( datname, get_data_filename() );
  strcpy( rootname, get_root_filename() );
  strcpy( rootname_loop, get_loop_filename() );
  strcpy( nsbpathname, get_nsb_directory() );
  strcpy( nsbpath_outer, get_nsb_directory_outer() );

  // get different parameters of the simulation

  addElecNoise = add_elec_noise(&FADC_noise_inner, &FADC_noise_outer, &DIGITAL_noise, &Trigger_noise);
  simulateNSB = get_nsb( &meanNSB, &nphe2NSB );  
  missPointing = get_misspointing(&missP_x,&missP_y);
  Spotsigma = get_sigma_xy_cm_spot(&Spot_x,&Spot_y);
  pmt_jitter = get_pmt_jitter_ns();


  // get selections on the parameters
  
  Select_Energy = get_select_energy( &Select_Energy_le, &Select_Energy_ue);
 
  //Parameters for starfield rotation

    Starfield_rotate = get_starfield_rotate();
 
  // log filenames information
  for(Int_t ict=0;ict<ct_Number;ict++){

    log(SIGNATURE,"\t%s : %i\n\t%20s:\t%s\n",
        "------- TELESCOPE ",ict+1,
        "Geometry ", GeometryName[ict]);
    strcpy( qe_filename, get_qe_filename(ict));

    printf("Telescope coordinates (cm): %.1f %.1f %.1f\n", CTx[ict], CTy[ict], CTz[ict]);

    // Look to factor for NSB respect to emi_coat PMTs
    if(strstr(qe_filename,"qe-emi-coat.RFL.dat") != 0)
      factorqe_NSB[ict]=EMICOAT_NSB;
    else if(strstr(qe_filename,"qe-emi.RFL.dat") != 0)
      factorqe_NSB[ict]=EMI_NSB;
    else if(strstr(qe_filename, "qe-intevac_hpd.RFL.dat") != 0)
      factorqe_NSB[ict]=HPD_INTEVAC_NSB;
    else if(strstr(qe_filename, "qe-hamamatsu_hpd.RFL.dat") != 0)
      factorqe_NSB[ict]=HPD_HAMAMATSU_NSB;
    else{
      log(SIGNATURE,"\n\nWARNING!! : the factor for the diffuse NSB for this QE file (%s) is not known to the camera simulation. The number of phe(ns stated in the inputcard will be used.\n\n",qe_filename); 
      factorqe_NSB[ict]=1.0;
    }
    log(SIGNATURE,
        "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s\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_CT[ict], 
        "Stars", starfieldname,
        "NSB database","(inner pixels)", nsbpathname,
        "(outer pixels)", nsbpath_outer,
        "QE", qe_filename, 
        "Data", datname,
        "ROOT",  rootname 
        );
    
  // log parameters information
  
  log(SIGNATURE,
      "%s:\n\t%20s: %f: %s\n\t%20s: %f\n\t%20s: %f\n",
      "Parameters",
      "NSB (phes/ ns 0.1*0.1 deg^2 239 m^2)", meanNSB*factorqe_NSB[ict], ONoff(simulateNSB),
      "Pedestals = ", get_FADC_pedestal(),
      "High to Low gain = ", FADC_high2low
      );
  
    // 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[ict][0],
          "Multiplicity",Trigger_multiplicity[ict],
          "Topology",Trigger_topology[ict]);
    }   
    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[ict],
          "Topology",Trigger_topology[ict]);
    }
    log(SIGNATURE,"\t%s\n",
        "END TELESCOPE --------");
  }
  // log flags information
  
  log(SIGNATURE,
      "%s:\n\t%23s: %s\n\t%23s: %s\n\t%23s: %s\n\t%23s: %s\n\t\t %s %3.2f,\n\t\t %s %3.2f (%s), %3.2f (%s) \n\t\t + %3.2f (%s)\n\n",
      "Flags",
      "Data_From_STDIN",   ONoff(Data_From_STDIN),  
      "Write_All_Event_Headers",  ONoff(Write_All_Event_Headers), 
      "Rotate Starfield",  ONoff(Starfield_rotate),
      "Electronic Noise",  ONoff(addElecNoise), 
      "Trigger noise (amplitude relative to that of single phe): ", 
      Trigger_noise, 
      "FADC noise (ADC counts): ", FADC_noise_inner, "inner p.",
      FADC_noise_outer, "outer p.", 
      DIGITAL_noise, "added digital noise at FADC"
      );

  if (! apply_gain_fluctuations())
    {
      log(SIGNATURE, "PMT gain fluctuations for the signal phes have been switched OFF");
      log(SIGNATURE, "in the input card.\n\n"); 
    }
  if (! apply_noise_gain_fluctuations())
    {
      log(SIGNATURE, "PMT gain fluctuations for the NSB phes have been switched OFF\n");
      log(SIGNATURE, "in the input card.\n\n"); 
    }

  if ( apply_gain_fluctuations() && !apply_noise_gain_fluctuations())
    {
      log(SIGNATURE, "Warning: you switched off PMT gain fluctuations for the NSB\n");
      log(SIGNATURE, "photoelectrons, but not for the signal photoelectrons.\n\n");
    }
  else if ( !apply_gain_fluctuations() && apply_noise_gain_fluctuations())
    {
      log(SIGNATURE, "Warning: you switched off PMT gain fluctuations for the signal\n");
      log(SIGNATURE, "photoelectrons, but not for the NSB photoelectrons.\n");
    }


  // log flags information

  log(SIGNATURE,
      "%s:\n\t%20s: %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 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)+1]; 
  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;
  }    

  // 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 (int i=0; i<nSkip; ++i)
      log(SIGNATURE, "\tshower # %d\n", Skip[i]);
  }

  
  // read parameters from the ct.def file

  //read_ct_file();

  Int_t Lev0, Lev1; 
  Int_t Lev0MT[MAX_NUMBER_OF_CTS], Lev1MT[MAX_NUMBER_OF_CTS]; 
  
  fadcValues    =  new TArrayC(FADC_slices_written);
  fadcValuesLow =  new TArrayC(FADC_slices_written);

  // number of pixels for parameters
    
  // Switched off writing TObject

  MArray::Class()->IgnoreTObjectStreamer();
  MParContainer::Class()->IgnoreTObjectStreamer();

  // initialise ROOT

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

  // initialise instance of Trigger and FADC classes

  MTrigger **Trigger_CT;
  Trigger_CT = new MTrigger *[ct_Number];

  for (int i=0; i<ct_Number;i++){
    Trigger_CT[i] = new MTrigger(TriggerPixels[i],
                    ((MGeomCam*)(camgeom.UncheckedAt(i))),
                    Trigger_gate_length, 
                    Trigger_overlaping_time,
                    Trigger_response_ampl, 
                    Trigger_response_fwhm, i);  //@< A instance of the Class MTrigger 
    Trigger_CT[i]->SetSeed(UInt_t(i+get_seeds(0)));

    if ( ! apply_gain_fluctuations())
      Trigger_CT[i]->SetGainFluctuations(kFALSE);
  }

  for (int ict=0; ict<ct_Number;ict++){
  // Generate database for trigger electronic noise
    if (addElecNoise)    
      Trigger_CT[ict]->SetElecNoise(Trigger_noise) ;

  // Set Right Discriminator threshold, taking into account trigger pixels
    Trigger_CT[ict]->CheckThreshold(&qThreshold[ict][0],GeometryCamera[ict]);
  }
  // Set flag in pixel 0 (not used for trigger) that indicates if secure pixel
  // is active: secureDiskThres*10000+riseDiskThres

  for (int ict=0; ict<ct_Number;ict++){
    if(riseDiskThres>0.0)
      qThreshold[ict][0]=((UInt_t)secureDiskThres*100)*100+riseDiskThres;
    else
      qThreshold[ict][0]=0.0;
  } 
  //  Initialise McTrig information class if we want to save trigger informtion

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

  int numberBranches;

  if (!Trigger_Loop)
    numberBranches=ct_Number;
  else
    numberBranches=icontrigger;

  if (Write_McTrig){

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

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

  if (Write_McFADC){

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

  float *fadc_jitter = new float[ct_Number];
  for (int i=0; i<ct_Number;i++)
    fadc_jitter[i] = 0.;

  MFadc **Fadc_CT;
  Fadc_CT = new MFadc *[ct_Number];

  for (int i=0; i<ct_Number;i++)
    Fadc_CT[i] = 
      new  MFadc(((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels(),
                 FADC_shape,
                 FADC_response_integ,FADC_response_fwhm,
                 FADC_shape_out,
                 FADC_resp_integ_out,FADC_resp_fwhm_out,
                 get_trig_delay(),
                 FADC_slices_per_ns,
                 FADC_slices_written); //@< 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[ct_NPixels];

  for(int i=0;i<ct_NPixels;i++)
    input_pedestals[i]=get_FADC_pedestal();
  for (int ict=0; ict<ct_Number;ict++){

    Fadc_CT[ict]->SetPedestals(input_pedestals);
    Fadc_CT[ict]->SetHigh2LowGain(FADC_high2low);
    
    // Generate database for the Fadc electronic noise

    if (!addElecNoise)
        continue;

    MGeomCam *camg = (MGeomCam*)(camgeom.UncheckedAt(ict));
    UInt_t n_inner_pixels;  // Number of inner(small) pixels
    for (n_inner_pixels = 0; n_inner_pixels < camg->GetNumPixels(); 
         n_inner_pixels ++)
    {
        if (camg->GetPixRatio(n_inner_pixels) < 1.)
            break;
    }
    Fadc_CT[ict]->SetElecNoise(FADC_noise_inner, FADC_noise_outer, 
                               n_inner_pixels);
    Fadc_CT[ict]->SetDigitalNoise(DIGITAL_noise);
  }

  // Prepare the raw data output

       // Header Tree
  MGeomCam **camdummy = new MGeomCam * [ct_Number];
  for (int ict=0; ict<ct_Number;ict++){
    camdummy[ict]= ((MGeomCam*)(camgeom.UncheckedAt(ict)));
  }
  MRawRunHeader *RunHeader= new MRawRunHeader();
  MMcRunHeader  *McRunHeader = new MMcRunHeader();
  MMcCorsikaRunHeader  *McCorsikaRunHeader = new MMcCorsikaRunHeader("","",ct_Number);


  for (int i = 0; i < ct_Number; i++)
    McCorsikaRunHeader->FillCT(CTx[i], CTy[i], CTz[i], -1., -1., -1., -1., i);


  MMcConfigRunHeader  **McConfigRunHeader = NULL;
  McConfigRunHeader = new MMcConfigRunHeader * [numberBranches];
  for (int i=0;i<numberBranches;i++) {
    McConfigRunHeader[i] = new MMcConfigRunHeader();
  }

       // Header branch

  MRawEvtHeader **EvtHeader = NULL;

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

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

  }

       // Data branch

  MRawEvtData **EvtData = NULL;    // Data branch

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

    for (int i=0;i<numberBranches;i++) {
      EvtData[i] = new MRawEvtData();
      EvtData[i]->InitRead(RunHeader);   //  We need the RunHeader to read
                                         //  number of pixels
    }
  }

  MMcEvt      **McEvt = NULL;
  MMcEvtBasic **McEvtBasic = NULL;

  if (Write_McEvt) 
    {
      McEvt = new MMcEvt *[ct_Number]; 
      McEvtBasic = new MMcEvtBasic *[ct_Number]; 
      for (int i = 0; i < ct_Number; i++)
        {
          McEvt[i] = new MMcEvt();
          McEvtBasic[i] = new MMcEvtBasic();
        }
    }
  //  
  // 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);

  HeaderTree.Branch("MMcCorsikaRunHeader.","MMcCorsikaRunHeader",
                    &McCorsikaRunHeader);

  if(!Trigger_Loop && Write_McTrig && ct_Number==1){
    
    HeaderTree.Branch("MMcTrigHeader.","MMcTrigHeader", 
                      &HeaderTrig[0]);
  }

  if (ct_Number==1){
    //    HeaderTree.Branch("MGeomCam.", "MGeomCam", &camdummy[0]);
    HeaderTree.Branch("MGeomCam.", GeometryName[0], &camdummy[0]);
    HeaderTree.Branch("MMcConfigRunHeader.","MMcConfigRunHeader",
                      &McConfigRunHeader[0]);
  }
  else{
    char branchname[256];
    for (int ict=0; ict<ct_Number;ict++){
      sprintf(help,"%i",ict+1);
      strcpy (branchname, "MGeomCam;");
      strcat (branchname, & help[0]);
      strcat (branchname, ".");
      HeaderTree.Branch(branchname, GeometryName[ict], 
                        &camdummy[ict]);
    }
    for(int i=0;i<numberBranches;i++){
      sprintf(help,"%i",i+1);
      strcpy (branchname, "MMcConfigRunHeader;");
      strcat (branchname, & help[0]);
      strcat (branchname, ".");
      HeaderTree.Branch(branchname,"MMcConfigRunHeader",
                        &McConfigRunHeader[i]);
    }
    
  }

  if ((Trigger_Loop || ct_Number>1) && Write_McTrig){
    ibr=0;
    for(char branchname[256];ibr<numberBranches;ibr++){
      
      sprintf(help,"%i",ibr+1);
      strcpy (branchname, "MMcTrigHeader;");
      strcat (branchname, & help[0]);
      strcat (branchname, ".");
      HeaderTree.Branch(branchname,"MMcTrigHeader", 
                        &HeaderTrig[ibr]);
    }
  }  

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

      HeaderTree.Branch(branchname,"MMcFadcHeader", 
                        &HeaderFadc[ibr]);
    }
  }  

  //  Fill branches for MRawRunHeader

  RunHeader->SetMagicNumber(MRawRunHeader::kMagicNumber);
  RunHeader->SetFormatVersion(5);
  RunHeader->SetSoftVersion((UShort_t) (VERSION*10));
  RunHeader->SetRunType(256);
  RunHeader->SetRunNumber(0);
  RunHeader->SetNumSamples(FADC_slices_written, FADC_slices_written);
  RunHeader->SetNumCrates(1);
  RunHeader->SetNumPixInCrate(ct_NPixels);

  //  Fill branches for MMcTrigHeader

  if(!Trigger_Loop && Write_McTrig && ct_Number==1){

    HeaderTrig[0]->SetTopology((Short_t) Trigger_topology[0]);
    HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity[0]);
    HeaderTrig[0]->SetThreshold(qThreshold[0]);
    HeaderTrig[0]->SetAmplitud(Trigger_response_ampl);
    HeaderTrig[0]->SetFwhm(Trigger_response_fwhm);
    HeaderTrig[0]->SetOverlap(Trigger_overlaping_time);
    HeaderTrig[0]->SetGate(Trigger_gate_length);
    HeaderTrig[0]->SetElecNoise(Trigger_noise);
    HeaderTrig[0]->SetGainFluctuations(Trigger_CT[0]->GetGainFluctuations());
    HeaderTrig[0]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
  }

  if(!Trigger_Loop && Write_McTrig && ct_Number>1){
    for(int i=0;i<ct_Number;i++){
      HeaderTrig[i]->SetTopology((Short_t) Trigger_topology[i]);
      HeaderTrig[i]->SetMultiplicity((Short_t) Trigger_multiplicity[i]);
      HeaderTrig[i]->SetThreshold(qThreshold[i]);
      HeaderTrig[i]->SetAmplitud(Trigger_response_ampl);
      HeaderTrig[i]->SetFwhm(Trigger_response_fwhm);
      HeaderTrig[i]->SetOverlap(Trigger_overlaping_time);
      HeaderTrig[i]->SetGate(Trigger_gate_length);
      HeaderTrig[i]->SetElecNoise(Trigger_noise);
      HeaderTrig[i]->SetGainFluctuations(Trigger_CT[i]->GetGainFluctuations());
      HeaderTrig[i]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
    }
  }  
  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<ct_NPixels;i++){
            fpixelthres[i]=
              ((Float_t)(fthrescount)>=qThreshold[0][i])?
              (Float_t)(fthrescount):qThreshold[0][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);
          HeaderTrig[iconcount]->SetElecNoise(Trigger_noise);
          HeaderTrig[iconcount]->SetGainFluctuations(Trigger_CT[0]->GetGainFluctuations());
          HeaderTrig[iconcount]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());
          iconcount++;
        }
      }   
    }
  }
  
  //  Fill branches for MMcFadcHeader
  Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);

  if(!Trigger_Loop && Write_McFADC && ct_Number==1){

      for(int k = 0; k < ct_NPixels; k++){
          if ( ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetPixRatio(k) < 1.)
              fadc_elecnoise[k]=FADC_noise_outer; // outer pixels
          else 
              fadc_elecnoise[k]=FADC_noise_inner; // inner pixels

          fadc_diginoise[k]=DIGITAL_noise;
      }

    HeaderFadc[0]->SetShape((Float_t)FADC_shape);
    HeaderFadc[0]->SetShapeOuter((Float_t)FADC_shape_out);
    HeaderFadc[0]->SetAmplitud(FADC_response_integ,
                               FADC_resp_integ_out);
    HeaderFadc[0]->SetFwhm(FADC_response_fwhm,FADC_resp_fwhm_out);
    HeaderFadc[0]->SetLow2High(FADC_high2low);
    HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],
                               ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
    HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
                                ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());

    //
    // Fill also the flag indicating whether the PMT gain fluctuations 
    // were simulated or not:
    //
    HeaderFadc[0]->SetGainFluctuations(Trigger_CT[0]->GetGainFluctuations());
    HeaderFadc[0]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());

  }

  if(!Trigger_Loop && Write_McFADC && ct_Number>1){
    for(int i=0;i<ct_Number;i++){

      for(int k = 0; k < ct_NPixels; k++){
          if ( ((MGeomCam*)(camgeom.UncheckedAt(i)))->GetPixRatio(k) < 1.)
              fadc_elecnoise[k]=FADC_noise_outer; // outer pixels
          else 
              fadc_elecnoise[k]=FADC_noise_inner; // inner pixels

          fadc_diginoise[k]=DIGITAL_noise;
      }

      Fadc_CT[i]->GetPedestals(&fadc_pedestals[0]);
      HeaderFadc[i]->SetShape((Float_t)FADC_shape);
      HeaderFadc[i]->SetShapeOuter((Float_t)FADC_shape_out);
      HeaderFadc[i]->SetAmplitud(FADC_response_integ,
                                 FADC_resp_integ_out);
      HeaderFadc[i]->SetFwhm(FADC_response_fwhm,FADC_resp_fwhm_out);
      HeaderFadc[i]->SetLow2High(FADC_high2low);
      HeaderFadc[i]->SetPedestal(&fadc_pedestals[0],((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
      HeaderFadc[i]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
                                  ((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());

      //
      // Fill also the flag indicating whether the PMT gain fluctuations 
      // were simulated or not:
      //
      HeaderFadc[i]->SetGainFluctuations(Trigger_CT[i]->GetGainFluctuations());
      HeaderFadc[i]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());

    }
  }
  
  if(Trigger_Loop && Write_McFADC){

    for(int k = 0; k < ct_NPixels; k++){
        if ( ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetPixRatio(k) < 1.)
            fadc_elecnoise[k]=FADC_noise_outer; // outer
        else 
            fadc_elecnoise[k]=FADC_noise_inner; // inner

        fadc_diginoise[k]=DIGITAL_noise;
    }

    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++){
          Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
          HeaderFadc[iconcount]->SetShape((Float_t)FADC_shape);
          HeaderFadc[iconcount]->SetShapeOuter((Float_t)FADC_shape_out);
          HeaderFadc[iconcount]->SetAmplitud(FADC_response_integ,
                                             FADC_resp_integ_out);
          HeaderFadc[iconcount]->SetFwhm(FADC_response_fwhm,
                                         FADC_resp_fwhm_out);
          HeaderFadc[iconcount]->SetLow2High(FADC_high2low);
          HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0], ct_NPixels);
          HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],
                                              &fadc_diginoise[0],ct_NPixels);

          HeaderFadc[iconcount]->SetGainFluctuations(Trigger_CT[0]->GetGainFluctuations());
          HeaderFadc[iconcount]->SetNoiseGainFluctuations((Bool_t)apply_noise_gain_fluctuations());

          iconcount++;
        }
      }   
    }
  }

  // UnSet flag in pixel 0 (not used for trigger) that indicates if secure pixel
  // is active once it has been stored

  for (int ict=0; ict<ct_Number;ict++){
      qThreshold[ict][0]=999999.99;
  } 

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

  if (Write_McEvt && ct_Number == 1)
    {
      EvtTree.Branch("MMcEvt.","MMcEvt", &McEvt[0]);
      EvtTree.Branch("MMcEvtBasic.","MMcEvtBasic", &McEvtBasic[0]);
    }

  if (Write_McEvt && ct_Number!=1)
    {
      ibr=0;
      for(char branchname[256]; ibr < numberBranches; ibr++)
        {
          sprintf(help,"%i",ibr+1);
          strcpy (branchname, "MMcEvt;");
          strcat (branchname, & help[0]);
          strcat (branchname, ".");
          EvtTree.Branch(branchname,"MMcEvt", 
                         &McEvt[ibr]);
        }
      for(char branchname[256]; ibr < numberBranches; ibr++)
        {
          sprintf(help,"%i",ibr+1);
          strcpy (branchname, "MMcEvtBasic;");
          strcat (branchname, & help[0]);
          strcat (branchname, ".");
          EvtTree.Branch(branchname,"MMcEvtBasic", 
                         &McEvtBasic[ibr]);
        }
    }

  if(!Trigger_Loop && ct_Number==1){
    
    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){
      ibr=0;
      for(char branchname[256];ibr<numberBranches;ibr++){
      
        sprintf(help,"%i",ibr+1);
        strcpy (branchname, "MMcTrig;");
        strcat (branchname, & help[0]);
        strcat (branchname, ".");
        EvtTree.Branch(branchname,"MMcTrig", 
                       &McTrig[ibr]);
      }
    }
  }  

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


  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;
      }
  }
  
  //  prepare the NSB simulation

  //  Instance of the Mlons class

  MLons lons(0, Trigger_response_ampl, Trigger_response_fwhm,
             FADC_shape, FADC_response_integ, FADC_response_fwhm,
             FADC_slices_per_ns, apply_noise_gain_fluctuations());

  lons.SetSeed(((UInt_t)(get_seeds(1)*get_seeds(1)*get_seeds(0))));

  lons.SetPath(nsbpathname);

  //  Instance of the Mlons class
  MLons lons_outer(0, Trigger_response_ampl, Trigger_response_fwhm,
                   FADC_shape_out,FADC_resp_integ_out,FADC_resp_fwhm_out,
                   FADC_slices_per_ns, apply_noise_gain_fluctuations());

  lons_outer.SetSeed(((UInt_t)(get_seeds(1)*get_seeds(0)*get_seeds(0))));

  lons_outer.SetPath(nsbpath_outer);

  if( simulateNSB){

    //
    // Calculate the non-diffuse NSB photoelectron rates
    //
    
    //
    // FIXME! --- star NSB different for each camera?
    // Then we will have to use mirror_frac[ict]
    //
    log(SIGNATURE,"Produce NSB rates from Star Field...\n");

    k = produce_nsbrates( starfieldname,
                          ((MGeomCam*)(camgeom.UncheckedAt(0))),
                          nsbrate_phepns,
                          0,
                          mirror_frac[0]);

    //
    // Call to "produce_nsbrates" above accounts ONLY for 
    // non-diffuse NSB (i.e. from stars). NOTE: produce_nsbrates already
    // accounts for the possibly different light collection efficiencies 
    // of inner and outer pixels, through a call (see function) to 
    // "produce_phes". The output array nsbrate_phepns contains only the 
    // rates due to individual stars in the FOV, and not diffuse NSB light!
    //

    if (k != 0){
      cout << "Error when reading starfield... \nExiting.\n";
      exit(1);
    }
  
    // calculate diffuse rate correcting for the pixel size and telescope

    float factorNSB_ct;
    for(int ict=0;ict<ct_Number;ict++){      
      
      // First we set the factor due to the mirror size with respect to the normal 
      // MAGIC mirror geometry.
      
      switch(GeometryCamera[ict]){
      case 1:
      case 2:
      case 3:
      case 4:
        factorNSB_ct=1.0;
        break;
      case 5:
      case 6:
      case 7:
        factorNSB_ct=1000.0/239.0;
        break;
      case 8:
      case 9:
      default:
        factorNSB_ct=1.0;
        break;  
      }

      factorNSB_ct = factorNSB_ct /
        ((*((MGeomCam*)(camgeom.UncheckedAt(ict)))).GetCameraDist()*1000*
         (*((MGeomCam*)(camgeom.UncheckedAt(ict)))).GetCameraDist()*1000*
         PIXEL_SIZE*PIXEL_SIZE); //  [mm^-2]

      for(UInt_t ui=0;
          ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); ui++){
        const Float_t size=
          (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD(); 
        // Returns distance [mm] between two paralel sides of pixel


        Float_t diffusensb = meanNSB*mirror_frac[ict];

        //
        // If pixel is outer pixel:
        //
        if ( ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetPixRatio(ui) < 1. )
          diffusensb *= WC_outer[1][90] / WC[1][90];
        //
        // FIXME! Correction above is for (possibly) different light collection efficiencies of 
        // inner and outer pixels. For the moment we assume the angular dependence of the light
        // collection is the same and hence use simply the ratio of efficiencies for light impinging
        // perpendicular to the camera plane (index 90 stands for 90 degrees)
        //

        diffnsb_phepns[ict][ui] = (Int_t(diffusensb*factorNSB_ct*100*size*size+0.5))/(100.0)*factorqe_NSB[ict];
      }
    }

    // calculate nsb rate including diffuse and starlight
    // we also include the extinction effect
    for(int ict=0;ict<ct_Number;ict++){      
      for(j=0;j<iNUMWAVEBANDS;j++){
        // calculate the effect of the atmospheric extinction (for stars!)
        
        zenfactor = pow(10., -0.4 * ext[j] );
        
        for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++)
          {
            nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + 
              zenfactor * nsbrate_phepns[ui][j];
            nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui];
          }
      }
    }
  }


  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  //
  //  Now 500 empty events with the condition in which the camera is run
  //  are simulated. In this way one gets an estimation of the 
  //  fluctuations of the pedestal of each FADC channel.
  //  This computation is done assuming any noise that affects 
  //  the FADC but there is no rotation of the Star Field (otherwise it
  //  should be done for each event). So it is valid if no starfield
  //  rotation is used.
  //
  //  Changed 20/03/2004, AM: now we no longer calculate the individual 
  //  FADC slice RMS. Due to correlations in the noise of neighboring 
  //  slides, it follows that RMS(sum_n_slices) != sqrt(n)*RMS(single_slice)
  //
  //  In the analysis in Mars, however, the RMS of the fluctuations of the 
  //  signal (resulting from the integration of n slices) is estimated as 
  //  sqrt(n) * MMcFadcHeader.fPedesSigmaHigh, where the latter value, 
  //  stored in the camera output, is calculated in the next lines of code.
  //  We have then made the following, as is being done also in real data: 
  //  calculate the RMS of the distribution of the sum of 14 slices, then
  //  the stored value is 
  //  MMcFadcHeader.fPedesSigmaHigh = RMS(sum_14_slices)/sqrt(14), and the
  //  same for the low gain. It can be seen that the RMS of the sum of n 
  //  slices, with n not too low (n>=6 or so), is more or less:
  //  RMS(sum_n_slices) ~ sqrt(n) * RMS(sum_14_slices)/sqrt(14)
  //
  //  The reason to sum 14 slices (and not 15) comes from the fact that in 
  //  real data there is a FADC clock noise affecting differently odd and 
  //  even FADC slices, so always an even number of them is added up so that
  //  this cancels out. So we do the calculation in the same way as in real 
  //  data.
  //
  //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

  Int_t empty_events = 500;
  Int_t n_slices = 14;

  for(int ict=0;ict<ct_Number;ict++){
    for(UInt_t ui=0;
        ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
        ui++){
      fadc_sigma[ui]=0.0;
      fadc_sigma_low[ui]=0.0;
    }
    for(int ie=0; ie < empty_events; ie++){
      Fadc_CT[ict]->Reset(); 
      if (addElecNoise){
        Fadc_CT[ict]->ElecNoise() ;
      }
      if(simulateNSB)
        {
          for(UInt_t ui=0;
              ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); 
              ui++)
            {
              if(nsb_phepns[ict][ui]>0.0)
                {
                  if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD() > 
                     (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD())
                    {

                      k = lons_outer.GetResponse(nsb_phepns[ict][ui],0.01,
                                                 & nsb_trigresp[0],
                                                 & nsb_fadcresp[0]);
                    }
                  else
                    {
                      k = lons.GetResponse(nsb_phepns[ict][ui],0.01,
                                           & nsb_trigresp[0],& nsb_fadcresp[0]);
                    }

                  if(k==0)
                    {
                      cout << "Exiting.\n";
                      exit(1);
                    }
                  Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
                }
            }
        }
      Fadc_CT[ict]->Pedestals();

      for(UInt_t ui=0;
          ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
          ui++)
        {
          //
          // We add up n_slices FADC slices (pedestal subtracted), then 
          // calculate the sigma_n-1 of this sum for the number of generated
          // noise events (=empty_events).
          //
          Float_t sumslices = Fadc_CT[ict]->AddNoiseInSlices(ui,1,n_slices);
          fadc_sigma[ui] += sumslices*sumslices;

          // Now the low gain:
          sumslices = Fadc_CT[ict]->AddNoiseInSlices(ui,0,n_slices);
          fadc_sigma_low[ui]+= sumslices*sumslices;
        }
    }

    for(UInt_t ui=0;
        ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
        ui++)
      {
        Float_t s_high = fadc_sigma[ui] / (Float_t)(empty_events-1) /
          (Float_t)(n_slices);
        fadc_sigma[ui] = sqrt(s_high);

        Float_t s_low = fadc_sigma_low[ui] / (Float_t)(empty_events-1) / 
          (Float_t)(n_slices);
        fadc_sigma_low[ui] = sqrt(s_low);
      }

    HeaderFadc[ict]->SetPedestalSigma(&fadc_sigma_low[0],&fadc_sigma[0],
                                      ((MGeomCam*)(camgeom.UncheckedAt(ict)))
                                      ->GetNumPixels());
  }


  if (is_calibration_run())
    {
      DoCalibration(Fadc_CT, Trigger_CT, camgeom, nsb_trigresp, nsb_fadcresp,
                    &lons, &lons_outer, nsb_phepns, addElecNoise, 
                    &EvtTree, EvtHeader, McEvt, EvtData);

      HeaderTree.Fill() ;

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

      cout << endl << "Calibration run finished!" << endl << endl;

      return (0);
    }

  //
  // Read the reflector file with the Cherenkov data
  //                    

  // select input file 

  if ( Data_From_STDIN )
    inputfile[0] = stdin;

  else
    {
      for(int ict = 0; ict < ct_Number; ict++)
        {
          log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] );
          inputfile[ict] = fopen( inname_CT[ict], "r" );
          
          if ( inputfile[ict] == NULL ) 
            error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] );
        }
    }
  
  // get signature, and check it

  for(int ict=0;ict<ct_Number;ict++)
    {
      if((reflector_file_version=check_reflector_file( inputfile[ict] ))==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" );
  strcpy( flag_new, "    \0" );

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

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

  // get flag

  for(int ict=0;ict<ct_Number;ict++)   
    fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
  

  // loop over the file

  still_in_loop = TRUE;

  // FIXME --- check if not eof for all input reflector files

  while (
         ((! Data_From_STDIN) && ( !feof(inputfile[0]) ))
         ||
         (Data_From_STDIN && still_in_loop)
         ) { 

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

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

      for(int ict=0;ict<ct_Number;ict++) {
        
        fread( flag_new, 4, 1, inputfile[ict] );

        if(!isA( flag_new, FLAG_START_OF_HEADER)){
          
          //  We break the main loop
          cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
          cout<<"         We break the main loop"<<endl;
          break;      
        }

        //      fread((char*)&mcrunh,(SIZE_OF_MCRUNHEADER)*sizeof(float), 1, inputfile[ict] );
        //  AM, changed reading method, 02/2004:

        mcrunh.read(inputfile[ict]);

      }
      
      nshow=0;

      for(int ict=0;ict<ct_Number;ict++){ 
        fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
      }
      while( isA( flag, FLAG_START_OF_EVENT   )){ // while there is a next event
        for(int ict=0;ict<ct_Number;ict++) {
          fread( flag_new, 4, 1, inputfile[ict] );
          if(!isA( flag_new, FLAG_EVENT_HEADER)){

            //  We break while events loop
            cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<" "<<ict<<endl;
            cout<<"         We break the while events loop"<<endl;
            break;      
          }
        }
        

        //
        // Clear Trigger and Fadc
        //
        for(int ict=0;ict<ct_Number;ict++) {
          Trigger_CT[ict]->Reset() ; 
          Trigger_CT[ict]->ClearFirst();
          Trigger_CT[ict]->ClearZero();
          Fadc_CT[ict]->Reset() ; 
        }
        
        ++nshow;
        if((nshow+ntshow+1)%100 == 1)
          log(SIGNATURE, "Event %d(+%d)\n", nshow, ntshow);
        
        // get MCEventHeader
        if (reflector_file_version<6){
          for(int ict=0;ict<ct_Number;ict++)
            {
              mcevth[ict].read(inputfile[ict]);
            }

        }
        else{
          for(int ict=0;ict<ct_Number;ict++)
            {
              mcevth_2[ict].read(inputfile[ict]);
            }
        }

        //
        // AM March 2004 simplified impact parameter calculation,
        // and allowed for correct estimate also for telescopes not
        // placed at Corsika's origin of coordinates (0.,0.).
        // FIXME: telescope coordinates still set by the user in the
        // input card, since they are not available in the reflector 
        // file!
        //
        // Calculate impact parameter as distance between telescope 
        // location and shower axis. In the previous implementation,
        // the definition was the distance between telescope axis and 
        // shower axis, but this has less physical meaning in general.
        // Light yield depends above all on distance to shower axis!
        // Of course for shower axis paralel to telescope the old and
        // the new definitions of impact parameter agree.
        // 


        // read the direction of the incoming shower
        // It is done only for one telescope since it is suposed
        // to be the same shower for all of them
        if (reflector_file_version<6){
          thetashw = mcevth[0].get_theta();
          phishw = mcevth[0].get_phi();
        }
        else{
          thetashw = mcevth_2[0].get_theta();
          phishw = mcevth_2[0].get_phi();
        }
        Zenith=thetashw; Azimutal=phishw;
        
        // calculate vector for shower
        
        l1 = sin(thetashw)*cos(phishw);
        m1 = sin(thetashw)*sin(phishw);
        n1 = cos(thetashw);

        //
        // Note, A.M.:
        // Attention! "core vector from mcevth*.get_core 
        // is the vector from the core to the telescope!
        //
        Float_t core2ct_x;
        Float_t core2ct_y;

        if (reflector_file_version<6)
          mcevth[0].get_core(&core2ct_x, &core2ct_y);   
        else
          mcevth_2[0].get_core(&core2ct_x, &core2ct_y); 

        //
        // Then true core position in Corsika's system is:
        //

        coreX = CTx[0] - core2ct_x;
        coreY = CTy[0] - core2ct_y;

        // 
        // FIXME: This may not work fine for CTs at z != 0 !!
        //
        for(int ict=0;ict<ct_Number;ict++)
          impactD[ict] = dist_r_P( CTx[ict], CTy[ict], CTz[ict], 
                                   l1, m1, n1, coreX, coreY, 0. );


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

        inumphe=0;

        for(int ict=0;ict<ct_Number;ict++){
         
          // Read first and last time and put inumphe_CT[0] to 0

          if (reflector_file_version<6)
            mcevth[ict].get_times(&arrtmin_ns,&arrtmax_ns);
          else
            mcevth_2[ict].get_times(&arrtmin_ns,&arrtmax_ns);
          
          inumphe_CT[ict]=0;


          // Obtain the FADC jitter of 1 FADC slice. This is a time to be added to the 
          // time of all photons in an event, before digitalization of the signal. It is
          // therefore the same time shift for all pixels in a CT. 

          fadc_jitter[ict] = 
            (1./Fadc_CT[ict]->GetFadcSlicesPerNanosec()) * RandomNumber; //ns

          // read the photons and produce the photoelectrons

          k = produce_phes( inputfile[ict],
                            ((MGeomCam*)(camgeom.UncheckedAt(ict))),
                            WAVEBANDBOUND1,
                            WAVEBANDBOUND6,
                            Trigger_CT[ict],   // will be changed by the function!
                            Fadc_CT[ict],      // will be changed by the function!
                            &inumphe_CT[ict], // important for later: the size of photoe[]
                            fnpix,    // will be changed by the function!
                            &ncph[ict],    // will be changed by the function!
                            &arrtmin_ns, // will be changed by the function!
                            &arrtmax_ns, // will be changed by the function!
                            ict,
                            mirror_frac[ict],
                            fadc_jitter[ict]);

          inumphe = (inumphe < inumphe_CT[ict])? inumphe_CT[ict] : inumphe;
          // inumphe will be the maximum number of phe in any of the telescopes


          if( k != 0 ){ // non-zero return value means error
            cout << "Exiting.\n";
            exit(1);
          }
        }
        
        // NSB simulation

        if(simulateNSB && inumphe >= nphe2NSB)
          {

            if(Starfield_rotate){
            
              // Introduction rho angle
            
              zenith = thetashw;
              azimutal = phishw;
              C1 = 0.48 * sin(zenith) - 0.87 * cos(zenith) * cos(azimutal);
              C3 = (0.87 * cos(zenith) - 0.48 * sin(zenith) * cos(azimutal));
              C2 = sqrt( sin(zenith) * sin(zenith) * sin(azimutal) * sin(azimutal) + 
                         C3 * C3 );
              rho = acos( C1/C2 );
           
              if ( sin(azimutal) < 0)
                rho = 2 * 3.14159 - rho;
              else
                rho = rho;

              rho = rho*180/3.14159;

              // Rotation of the NSB
              // FIXME --- We should rotate for all cameras. Is it always the same rho?         
              for(int ict=0;ict<ct_Number;ict++)
                k = size_rotated(&nsb_phepns_rotated[ict][0],
                                 nsb_phepns[ict],
                                 rho);
            }

          //  Fill trigger and fadc response in the trigger class from the database
          for(int ict=0;ict<ct_Number;ict++)
            {
              for(UInt_t ui=0;
                  ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
                  ui++)
                {
                  if(nsb_phepns_rotated[ict][ui]>0.0)
                    {
                      if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD() > 
                         (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD())
                        {
                          k=lons_outer.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
                                                   & nsb_trigresp[0],
                                                   & nsb_fadcresp[0]);
                        }
                      else
                        {
                          k=lons.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
                                             & nsb_trigresp[0],& nsb_fadcresp[0]);
                        }
                      if(k==0)
                        {
                          cout << "Exiting.\n";
                          exit(1);
                        }
                      Trigger_CT[ict]->AddNSB(ui,nsb_trigresp);
                      Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
                    }
                }
            }
          
          }// end if(simulateNSB && inumphe_CT[0] >= nphe2NSB) ...


        for(int ict=0;ict<ct_Number;ict++)
          {
            inumphensb[ict]=0;

            for (UInt_t ui=0;
                 ui < ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
                 ui++)
              inumphensb[ict]+=nsb_phepns[ict][ui]*TOTAL_TRIGGER_TIME;

            ntcph[ict]+=ncph[ict];
            if ((nshow+ntshow+1)%100 == 1){
              log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
                  ncph[ict], ntcph[ict]);

              cout << "Total number of phes in CT "<<ict<<": " 
                   << inumphe_CT[ict]<<" (+ ";
              cout<<inumphensb[ict]<<" mean expected number from NSB)"<<endl;
            }
          }

        // skip it ?

        int i;  
        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;
        }
        
        //++++++++++++++++++++++++++++++++++++++++++++++++++
        // 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

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

        //   now the noise of the electronic 
        //   (preamps, optical transmission,..)  is introduced. 
        //   This is done inside the class MTrigger by the method ElecNoise. 
        //   

        for(int ict=0;ict<ct_Number;ict++) 
          {
            if (addElecNoise && inumphe >= nphe2NSB)
              { 
                Trigger_CT[ict]->ElecNoise(Trigger_noise);
                Fadc_CT[ict]->ElecNoise();
              }
          }
        
        //  now a shift in the fadc signal due to the pedestals is 
        //  introduced
        //  This is done inside the class MFadc by the method Pedestals

        for(int ict=0;ict<ct_Number;ict++)
          Fadc_CT[ict]->Pedestals();

        
        //   We study several trigger conditons
        if(Trigger_Loop)
          {
            // Set to zero the flag to know if some conditon has triggered
            btrigger=0;
            flagstoring = 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 (int i=0;i<ct_NPixels;i++)
                  {
                    fpixelthres[i] =
                      ((Float_t)(fthrescount)>=qThreshold[0][i])?
                      (Float_t)(fthrescount):qThreshold[0][i];

                    // Rise the discrimnator threshold to avoid huge rates

                    if(riseDiskThres>0.0 && simulateNSB && inumphe >= nphe2NSB)
                      for(int ii=0;ii<ct_NPixels;ii++)
                        if( nsb_phepns_rotated[0][ii]>riseDiskThres)
                          fpixelthres[ii]=secureDiskThres;
                  }
                Trigger_CT[0]->SetThreshold(fpixelthres);
            
                Trigger_CT[0]->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
                Lev0=0;
                Lev1=0;

                //  loop over multiplicity of trigger configuration
                for (imulticount = Trigger_loop_lmult;
                     imulticount <= Trigger_loop_umult;
                     imulticount++)
                  {
                    Trigger_CT[0]->SetMultiplicity(imulticount);
                    Trigger_CT[0]->ClearZero();

                    Lev0=(Short_t) Trigger_CT[0]->ZeroLevel();
                    if (Lev0>0 || Write_All_Event_Headers || btrigger)
                      {
                
                        //  loop over topologies
                        for(itopocount=Trigger_loop_ltop;
                            itopocount<=Trigger_loop_utop;
                            itopocount++)
                          {
                            Lev1=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_CT[0]->SetTopology(1);
                            else
                              {
                                // We should be careful that topologies are sort from 
                                // the less to the more restrictive one. 
                                Trigger_CT[0]->SetTopology(isorttopo[itopocount]);
                              }
                            Trigger_CT[0]->ClearFirst();
                  
                            //
                            //   Start the First Level Trigger simulation
                            //
                            if(Lev0!=0)
                              Lev1=Trigger_CT[0]->FirstLevel();
                            if(Lev1>0) {
                              btrigger= 1;
                              ntriggerloop[ithrescount]
                                [imulticount-Trigger_loop_lmult]
                                [itopocount-Trigger_loop_ltop]++;
                            }
                  
                            Lev0=1;
                            Int_t NumImages = Lev1;
                            if(Lev1==0 && (Write_All_Event_Headers || btrigger))
                              {
                                btrigger= 1;
                                NumImages=1;
                                Lev0=0;
                              }

                            for (Int_t ii=0;ii<NumImages;ii++)
                              {
                                if (Write_McTrig)
                                  {
                                    McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
                                    McTrig[iconcount]->
                                      SetTime(Trigger_CT[0]->GetFirstLevelTime(ii),ii+1);
                                    Trigger_CT[0]->GetMapDiskriminator(trigger_map);
                                    McTrig[iconcount]->SetMapPixels(trigger_map,ii);
                                  }
                                //
                                //  fill inside the class fadc the member output
                                //

                                Fadc_CT[0]->TriggeredFadc(Trigger_CT[0]->
                                                          GetFirstLevelTime(ii));
            
                                if( Write_RawEvt )
                                  {
                                    //
                                    //  Fill the header of this event 
                                    //
                      
                                    EvtHeader[iconcount]->
                                      FillHeader( (UInt_t) (ntshow + nshow),0);
                                    EvtHeader[iconcount]->
                                      SetTriggerPattern((UInt_t)MTriggerPattern::kTriggerLvl1);
                                    //   fill pixel information
                                    if (Lev1 /*|| Write_All_Event_Headers AMTEST*/){
                                      if (addElecNoise) Fadc_CT[0]->DigitalNoise();
                                      for(UInt_t i=0;
                                          i<((MGeomCam*)(camgeom.UncheckedAt(0)))
                                            ->GetNumPixels();i++){
                                        //
                                        // AM 15 01 2004: commented out "continue" 
                                        // statement, so that also pixels with no 
                                        // C-photons will be written to the output 
                                        // in case the camera is run with no noise.
                                        // if(!Fadc_CT[0]->IsPixelUsed(i)) continue;

                                        for (j=0;j<FADC_slices_written;j++){
                                          fadcValues->AddAt(Fadc_CT[0]->
                                                            GetFadcSignal(i,j),j);
                                          fadcValuesLow->AddAt(Fadc_CT[0]->
                                                               GetFadcLowGainSignal(i,j),j);
                                        }
                                        EvtData[iconcount]->AddPixel(i,fadcValues,0);
                                        EvtData[iconcount]->AddPixel(i,fadcValuesLow,kTRUE);
                                      }
                                    }
                                  }
                              }
                            //
                            // Increase counter of analised trigger conditions
                            //
                            iconcount++;
                          }
                      }
                    else{
                      break;
                    }
                  }
                if (!btrigger) break;
              }
            if (btrigger){

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

              if(!flagstoring)
                nstoredevents++;
              flagstoring = 1;
              
              if (Write_McEvt) {
                Float_t ftime, ltime;
                if (reflector_file_version<6){
                  mcevth[0].get_times(&ftime, &ltime);

                  McEvtBasic[0]->Fill((MMcEvt::ParticleId_t) mcevth[0].get_primary(),
                                      mcevth[0].get_energy(), impactD[0], 
                                      phiCT[0], thetaCT[0]);

                  McEvt[0]->Fill( 0, 
                                  (UShort_t) mcevth[0].get_primary() , 
                                  mcevth[0].get_energy(),
                                  -1.0,
                                  -1.0,
                                  -1.0, 
                                  mcevth[0].get_theta(), 
                                  mcevth[0].get_phi(), 
                                  mcevth[0].get_core(),
                                  coreX,
                                  coreY,
                                  impactD[0],
                                  phiCT[0],
                                  thetaCT[0],
                                  ftime,
                                  ltime,
                                  0,
                                  0,
                                  0,
                                  0,
                                  0,
                                  0,
                                  0,
                                  (UInt_t)mcevth[0].get_CORSIKA(), 
                                  (UInt_t)mcevth[0].get_AtmAbs(), 
                                  (UInt_t)(mcevth[0].get_MirrAbs()+
                                           mcevth[0].get_OutOfMirr()+
                                           mcevth[0].get_BlackSpot()), 
                                  (UInt_t) ncph[0], 
                                  (UInt_t) inumphe_CT[0],
                                  (UInt_t) inumphensb[0]+inumphe_CT[0],
                                  -1.0,
                                  -1.0,
                                  -1.0,
                                  fadc_jitter[0]);
                }
                else{
                  Float_t Nmax, t0, tmax, a, b, c, chi2;
                  mcevth_2[0].get_times(&ftime, &ltime);
                  chi2=mcevth_2[0].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);

                  McEvtBasic[0]->Fill((MMcEvt::ParticleId_t) mcevth_2[0].get_primary(),
                                      mcevth_2[0].get_energy(), impactD[0],
                                      mcevth_2[0].get_phi_CT(),
                                      mcevth_2[0].get_theta_CT());

                  McEvt[0]->Fill((UInt_t) mcevth_2[0].get_evt_number(),
                                 (UShort_t) mcevth_2[0].get_primary() , 
                                 mcevth_2[0].get_energy(),
                                 mcevth_2[0].get_thick0(),
                                 mcevth_2[0].get_first_target(),
                                 mcevth_2[0].get_z_first_int(),
                                 mcevth_2[0].get_theta(), 
                                 mcevth_2[0].get_phi(), 
                                 mcevth_2[0].get_core(),
                                 coreX,
                                 coreY,
                                 impactD[0],
                                 mcevth_2[0].get_phi_CT(),
                                 mcevth_2[0].get_theta_CT(),
                                 ftime,
                                 ltime,
                                 Nmax,
                                 t0,
                                 tmax,
                                 a,
                                 b,
                                 c,
                                 chi2,
                                 (UInt_t)mcevth_2[0].get_CORSIKA(), 
                                 (UInt_t)mcevth_2[0].get_AtmAbs(), 
                                 (UInt_t)(mcevth_2[0].get_MirrAbs()+
                                          mcevth_2[0].get_OutOfMirr()+
                                          mcevth_2[0].get_BlackSpot()), 
                                 (UInt_t) ncph[0], 
                                 (UInt_t) inumphe_CT[0],
                                 (UInt_t) inumphensb[0]+inumphe_CT[0],
                                 mcevth_2[0].get_ElecFraction(),
                                 mcevth_2[0].get_MuonFraction(),
                                 mcevth_2[0].get_OtherFraction(),
                                 fadc_jitter[0]);
                }
              }
              //  Fill the Tree with the current leaves of each branch
              i=EvtTree.Fill() ;

              //  Clear the branches
              if(Write_McTrig){
                for(int i=0;i<numberBranches;i++){
                  McTrig[i]->Clear() ;
                }
              }
              if( Write_RawEvt ){
                for(int i=0;i<numberBranches;i++){
                  EvtHeader[i]->Clear() ;
                  EvtData[i]->ResetPixels (0, 0);
                }
              }
              if (Write_McEvt)
                {
                  McEvt[0]->Clear();
                  McEvtBasic[0]->Clear();
                }
            }
          }

        //  We study a single trigger condition
        else {
          
          // Set to zero the flag to know if some conditon has triggered
          btrigger=0;
          flagstoring = 0;
          
          for(int ict = 0; ict < ct_Number; ict++)
            {

              //  Setting trigger conditions
              Trigger_CT[ict]->SetMultiplicity(Trigger_multiplicity[ict]);
              Trigger_CT[ict]->SetTopology(Trigger_topology[ict]);
              for (int i=0;i<ct_NPixels;i++)
                fpixelthres[i]=qThreshold[ict][i];

              // Rise the discrimnator threshold to avoid huge rates
              if(riseDiskThres>0.0 && simulateNSB && inumphe >= nphe2NSB)
                for(int ii=0;ii<ct_NPixels;ii++)
                  {
                    if( nsb_phepns_rotated[ict][ii]>riseDiskThres)
                      fpixelthres[ii]=secureDiskThres;
                  }
          
              Trigger_CT[ict]->SetThreshold(fpixelthres);
                                 
              Trigger_CT[ict]->Diskriminate() ; 

              //
              //   Look if in all the signals in the trigger signal branch
              //   is a possible Trigger. Therefore we have to discriminate all
              //   the simulated analog signals (Method Diskriminate in class
              //   MTrigger). We look simultaneously for the moments at which
              //   there are more than TRIGGER_MULTI pixels above the 
              //   CHANNEL_THRESHOLD. 
              //
            
            Lev0MT[ict] = (Short_t) Trigger_CT[ict]->ZeroLevel() ; 

            Lev1MT[ict] = 0 ; 
          
            //
            //   Start the First Level Trigger simulation
            //
          
            if ( Lev0MT[ict] > 0 /*|| Write_All_Event_Headers*/)
              Lev1MT[ict]= Trigger_CT[ict]->FirstLevel();

            if (Lev1MT[ict]>0)
              ++ntrigger[ict];

          }

          Int_t NumImages = 0;
          Int_t CT_triggered=0;
          for(int ict=0;ict<ct_Number;ict++)
            {
              if(NumImages==0 && Lev1MT[ict]>0)
                CT_triggered=ict;
              NumImages = (NumImages>=Lev1MT[ict]) ? NumImages : 1;
              Lev0MT[ict]=1;

              if (Lev1MT[ict]==0 && Write_All_Event_Headers)
                {
                  NumImages=1;
                  Lev0MT[ict]=0;
                }
            }

          for(int ict=0;ict<ct_Number;ict++){
            for(Int_t ii=0;ii<NumImages;ii++){

              btrigger=1;
              
              //  Loop over different level one triggers
              
              //
              //   fill inside class fadc the member output
              //
              if(Lev1MT[ict]>0)
                Fadc_CT[ict]->
                  TriggeredFadc(Trigger_CT[ict]->GetFirstLevelTime(ii));
              else
                Fadc_CT[ict]->
                  TriggeredFadc(Trigger_CT[CT_triggered]->GetFirstLevelTime(ii));

              if(!flagstoring)
                nstoredevents++;
              flagstoring = 1;

              if (Write_McTrig){
                McTrig[ict]->SetFirstLevel (Lev1MT[ict]);
                McTrig[ict]
                  ->SetTime(Trigger_CT[ict]->GetFirstLevelTime(ii),ii+1);
                Trigger_CT[ict]->GetMapDiskriminator(trigger_map);
                McTrig[ict]->SetMapPixels(trigger_map,ii);
              }
              
              //  Fill Evt information
              
              if (Write_RawEvt){
                
                //
                //  Fill the header of this event 
                //

                EvtHeader[ict]
                  ->FillHeader ( (UInt_t) (ntshow + nshow) , 0 ) ; 
                EvtHeader[ict]->SetTriggerPattern((UInt_t)MTriggerPattern::kTriggerLvl1);

                // Fill pixel information
                // AM 17/2/2005: added condition on inumphe. Noise is not generated for 
                // events with less than phe2NSB photoelectrons, and then it is better not
                // to write the images of those events to the output, so that we avoid them
                // to be processed (with no noise) later in the chain.
                
                if (Lev1MT[ict] && inumphe >= nphe2NSB)
                  {
                    if (addElecNoise) 
                      Fadc_CT[ict]->DigitalNoise();
                    for(UInt_t i=0;
                        i<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
                        i++){
//
// AM 15 01 2004 Commented out "continue" statement, so that also pixels
// containing no C-photons will be written to the output in case of running
// camera with no noise added to the signal.
//                  if(!Fadc_CT[ict]->IsPixelUsed(i)) continue;
//
                      for (j = 0; j < FADC_slices_written; j++)
                        {
                          fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(i,j),j);
                          fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(i,j),j);
                        }
                      EvtData[ict]->AddPixel(i,fadcValues,0);
                      EvtData[ict]->AddPixel(i,fadcValuesLow,kTRUE);
                    }
                  }
              } 
            }
            //
            //    if a first level trigger occurred, then 
            //       1. do some other stuff (not implemented) 
            //       2. start the gui tool
            
            if(FADC_Scan){
              if ( Lev0MT[ict] > 0 ) {
                  Fadc_CT[ict]->ShowSignal( McEvt[ict], (Float_t) 60. ) ; 
              }
            }
            
            if(Trigger_Scan){
              if ( Lev0MT[ict] > 0 ) {
                Trigger_CT[ict]->ShowSignal(McEvt[ict]) ;
              }
            }
              
          } // end CT loop

          // If there is trigger in some telescope or we store all showers
          if(btrigger){  
            if (Write_McEvt){
              //
              //   fill the MMcEvt with all information  
              //
            
              for (int ict=0;ict<ct_Number;ict++){
                Float_t ftime, ltime;
                if (reflector_file_version<6){
                  mcevth[ict].get_times(&ftime, &ltime);

                  McEvtBasic[ict]->Fill((MMcEvt::ParticleId_t) mcevth[0].get_primary(),
                                        mcevth[ict].get_energy(), impactD[ict],
                                        phiCT[ict], thetaCT[ict]);

                  McEvt[ict]->Fill( 0, 
                                  (UShort_t) mcevth[0].get_primary() , 
                                  mcevth[ict].get_energy(),
                                  -1.0,
                                  -1.0,
                                  -1.0, 
                                  mcevth[ict].get_theta(), 
                                  mcevth[ict].get_phi(), 
                                  mcevth[ict].get_core(),
                                  coreX,
                                  coreY,
                                  impactD[ict],
                                  phiCT[ict],
                                  thetaCT[ict],
                                  ftime,
                                  ltime,
                                  0,
                                  0,
                                  0,
                                  0,
                                  0,
                                  0,
                                  0,
                                  (UInt_t)mcevth[ict].get_CORSIKA(), 
                                  (UInt_t)mcevth[ict].get_AtmAbs(), 
                                  (UInt_t)(mcevth[ict].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()), 
                                  (UInt_t) ncph[ict], 
                                  (UInt_t) inumphe_CT[ict],
                                  (UInt_t) inumphensb[ict]+inumphe_CT[ict],
                                  -1.0,
                                  -1.0,
                                  -1.0,
                                   fadc_jitter[ict]);
                }
                else{
                  Float_t Nmax, t0, tmax, a, b, c, chi2;
                  mcevth_2[ict].get_times(&ftime, &ltime);
                  chi2=mcevth_2[ict].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);

                  McEvtBasic[ict]->Fill((MMcEvt::ParticleId_t)mcevth_2[ict].get_primary(), 
                                        mcevth_2[ict].get_energy(), impactD[ict],
                                        mcevth_2[ict].get_phi_CT(),
                                        mcevth_2[ict].get_theta_CT());
 

                  McEvt[ict]->Fill( (UInt_t) mcevth_2[ict].get_evt_number(),
                                  (UShort_t) mcevth_2[ict].get_primary() , 
                                  mcevth_2[ict].get_energy(),
                                  mcevth_2[ict].get_thick0(),
                                  mcevth_2[ict].get_first_target(),
                                  mcevth_2[ict].get_z_first_int(),
                                  mcevth_2[ict].get_theta(), 
                                  mcevth_2[ict].get_phi(), 
                                  mcevth_2[ict].get_core(),
                                  coreX,
                                  coreY,
                                  impactD[ict],
                                  mcevth_2[ict].get_phi_CT(),
                                  mcevth_2[ict].get_theta_CT(),
                                  ftime,
                                  ltime,
                                  Nmax,
                                  t0,
                                  tmax,
                                  a,
                                  b,
                                  c,
                                  chi2,
                                  (UInt_t)mcevth_2[ict].get_CORSIKA(), 
                                  (UInt_t)mcevth_2[ict].get_AtmAbs(), 
                                  (UInt_t) (mcevth_2[ict].get_MirrAbs()+mcevth_2[ict].get_OutOfMirr()+mcevth_2[ict].get_BlackSpot()), 
                                  (UInt_t) ncph[ict], 
                                  (UInt_t) inumphe_CT[ict],
                                  (UInt_t) inumphensb[ict]+inumphe_CT[ict],
                                  mcevth_2[ict].get_ElecFraction(),
                                  mcevth_2[ict].get_MuonFraction(),
                                  mcevth_2[ict].get_OtherFraction(),
                                  fadc_jitter[ict]);
                }
              }
            }
            //   We do not count photons out of the camera.     
              
            
            //
            //    write it out to the file outfile
            // 
            
            EvtTree.Fill() ; 
          }
            
          //    clear all
          for(int ict=0;ict<ct_Number;ict++){
            if (Write_RawEvt) EvtHeader[ict]->Clear() ;
            if (Write_RawEvt) EvtData[ict]->ResetPixels(0,0);
            if (Write_McTrig) McTrig[ict]->Clear() ;
            if (Write_McEvt) 
              {
                McEvt[ict]->Clear() ; 
                McEvtBasic[ict]->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<((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels(); ++i) {
          printf("%d (%d): ", i, npixneig[i]);
          for (j=0; j<npixneig[i]; ++i) 
            printf(" %d", pixneig[i][j]);
          printf("\n");
        }
        
#endif // __DEBUG__
                          

        // We search the maximum impact parameter fo the simualted showers
        maxpimpact=maxpimpact<impactD[0]?impactD[0]:maxpimpact;
        
        // look for the next event
        
        for(int ict=0;ict<ct_Number;ict++)
          fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );

      } // 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 );
        break;
      }
      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 );
        for(int ict=0;ict<ct_Number;ict++)
          fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
        
        if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
          log(SIGNATURE, "End of file . . .\n");
          still_in_loop  = FALSE;
          log(SIGNATURE, "Reading ASCII files at the end of the reflector file.  . .\n");
          for(int ict=0;ict<ct_Number;ict++){
            read_ascii(inputfile[ict], McConfigRunHeader[ict]);
            McConfigRunHeader[ict]->SetMissPointingX(missP_x);
            McConfigRunHeader[ict]->SetMissPointingY(missP_y);

            McConfigRunHeader[ict]->SetMirrorFraction(mirror_frac[ict]);

            if ( Spotsigma > 0.)
            {
                Float_t ref_spotsigma = McConfigRunHeader[ict]->GetPointSpread();
                Float_t newsigma = sqrt(ref_spotsigma * ref_spotsigma + 
                    Spot_x* Spot_x);
                McConfigRunHeader[ict]->SetPointSpreadX(newsigma);
                newsigma = sqrt(ref_spotsigma * ref_spotsigma + Spot_y* Spot_y);
                McConfigRunHeader[ict]->SetPointSpreadY(newsigma);
            }
          }
          if ((! Data_From_STDIN) && ( !feof(inputfile[0]) )){
            
            // we have concatenated input files.
            // get signature of the next part and check it.

            if((reflector_file_version=check_reflector_file( inputfile[0] ))==FALSE){
              log(SIGNATURE, "Next file is not recognised as a reflector file.\n");
              log(SIGNATURE, "Stopping ...\n");
              break;
            }
            
          }

        for(int ict=0;ict<ct_Number;ict++)        
          fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );

        } // 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 the McRunHeader

  Float_t heights[10];
  time_t ltime;
  Float_t ftime;
  Float_t rnum;
  Float_t  viewcone[2]={0,0};

  get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
  if (reflector_file_version<6){
    mcevth[0].get_theta_range(&shthetamin, &shthetamax);
    mcevth[0].get_phi_range(&shphimin,&shphimax);
    mcevth[0].get_theta_range(&shthetamin, &shthetamax);
    mcevth[0].get_phi_range(&shphimin,&shphimax);
    corsika=UInt_t(mcevth[0].get_VersionPGM()*1000);
    for (int i=0; i< 10;i++)
      heights[i]=mcevth[0].get_HeightLev (i);
    rnum=mcevth[0].get_RunNumber();
  }
  else{
    mcevth_2[0].get_theta_range(&shthetamin, &shthetamax);
    mcevth_2[0].get_phi_range(&shphimin,&shphimax);
    corsika=UInt_t(mcevth_2[0].get_VersionPGM()*1000);
    for (int i=0; i< 10;i++)
      heights[i]=mcevth_2[0].get_HeightLev (i);
    rnum=mcevth_2[0].get_RunNumber();
    mcevth_2[0].get_viewcone(&viewcone[0],&viewcone[1]);
  }

  if(!Trigger_Loop) icontrigger=0;
  time (&ltime);
  ftime = ((Float_t)ltime)/1000;

  if (reflector_file_version<6)
    McRunHeader->Fill(rnum,
                    (UInt_t) 0,
                    mcevth[0].get_DateRun(),
                    ftime,
                    icontrigger,
                    !Write_All_Event_Headers,
                    Write_McEvt,
                    Write_McTrig,
                    Write_McFADC,
                    Write_RawEvt,
                    addElecNoise,
                    ct_NPixels,
                    (UInt_t)ntshow,
                    (UInt_t)nstoredevents,
                    0,
                    sfRaH,
                    sfRaM,
                    sfRaS,
                    sfDeD,
                    sfDeM,
                    sfDeS,
                    meanNSB,
                    shthetamax,
                    shthetamin,
                    shphimax,
                    shphimin,
                    maxpimpact,
                    mcevth[0].get_CWaveLower(),
                    mcevth[0].get_CWaveUpper(),
                    mcevth[0].get_slope(),
                    1,
                    heights,
                    corsika,
                    (UInt_t)(reflector_file_version*100),
                    (UInt_t)(VERSION*100),
                    0);
  else
    McRunHeader->Fill(rnum,
                    (UInt_t) 0,
                    mcevth_2[0].get_DateRun(),
                    ftime,
                    icontrigger,
                    !Write_All_Event_Headers,
                    Write_McEvt,
                    Write_McTrig,
                    Write_McFADC,
                    Write_RawEvt,
                    addElecNoise,
                    ct_NPixels,
                    (UInt_t)ntshow,
                    (UInt_t)nstoredevents,
                    0,
                    sfRaH,
                    sfRaM,
                    sfRaS,
                    sfDeD,
                    sfDeM,
                    sfDeS,
                    meanNSB,
                    shthetamax,
                    shthetamin,
                    shphimax,
                    shphimin,
                    maxpimpact,
                    mcevth_2[0].get_CWaveLower(),
                    mcevth_2[0].get_CWaveUpper(),
                    mcevth_2[0].get_slope(),
                    1,
                    heights,
                    corsika,
                    (UInt_t)(reflector_file_version*100),
                    (UInt_t)(VERSION*100),
                    0);
  // Fill some missing values for MRawRunHeader

  RunHeader->SetRunNumber((UInt_t)rnum);
  RunHeader->SetNumEvents(nstoredevents);

  // Fill MMcCorsikaRunHeader

  Float_t constantC[50];
  mcrunh.get_constantC(&constantC[0]);
  Float_t constantCKA[40];
  mcrunh.get_constantCKA(&constantCKA[0]);
  Float_t constantCETA[5];
  mcrunh.get_constantCETA(&constantCETA[0]);
  Float_t constantCSTRBA[11];
  mcrunh.get_constantCSTRBA(&constantCSTRBA[0]);
  Float_t constantAATM[5];
  mcrunh.get_constantAATM(&constantAATM[0]);
  Float_t constantBATM[5];
  mcrunh.get_constantBATM(&constantBATM[0]);
  Float_t constantCATM[5];
  mcrunh.get_constantCATM(&constantCATM[0]);
  Float_t constantNFL[4];
  mcrunh.get_constantNFL(&constantNFL[0]);
  
  if(reflector_file_version>5)
    McCorsikaRunHeader->Fill(rnum,
                             mcrunh.get_date(),
                             corsika,
                             1,
                             heights,
                             mcevth_2[0].get_slope(),
                             mcrunh.get_ELow(),
                             mcrunh.get_EUpp(),
                             mcrunh.get_EGS4(),
                             mcrunh.get_NKG(),
                             mcrunh.get_Ecutoffh(),
                             mcrunh.get_Ecutoffm(),
                             mcrunh.get_Ecutoffe(),
                             mcrunh.get_Ecutoffg(),
                             constantC,
                             constantCKA,
                             constantCETA,
                             constantCSTRBA,
                             constantAATM,
                             constantBATM,
                             constantCATM,
                             constantNFL,
                             viewcone,
                             mcrunh.get_wobble(),
                             mcrunh.get_atmophere()
                             );
  
    //  Store qe for each PMT in output file
    TArrayF qe_pmt;
    TArrayF wav_pmt;

  for(int ict=0;ict<ct_Number;ict++){
    McConfigRunHeader[ict]->InitSizePMTs(ct_NPixels);
    for(int i=0; i<(Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();i++){
        McConfigRunHeader[ict]->SetPmtTimeJitter(pmt_jitter);
        McConfigRunHeader[ict]->AddPMT(i);
        MGeomPMT &pmt  = McConfigRunHeader[ict]->GetPMT(i);
        qe_pmt.Set(pointsQE[ict],QE[ict][i][1]);
        wav_pmt.Set(pointsQE[ict],QElambda);
        pmt.SetArraySize(pointsQE[ict]);
        pmt.SetPMTContent(i,wav_pmt,qe_pmt);
    }
    
    // Store Light Collection factors in the output file
    TArrayF theta_lc;
    TArrayF factor_lc;
    TArrayF factor_lc_outer;

    theta_lc.Set(pointsWC,WC[0]);
    factor_lc.Set(pointsWC,WC[1]);
    factor_lc_outer.Set(pointsWC,WC_outer[1]);

    McConfigRunHeader[ict]->SetLightCollection(theta_lc, factor_lc, 
                                               factor_lc_outer);
  
  }

  //  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
  
  if (Trigger_Loop){
    log( SIGNATURE, "%d event(s), with a total of %d C.photons\n", 
         ntshow, ntcph[0] );
    datafile<<ntshow<<" event(s), with a total of "<<ntcph[0]<<" C.photons"<<endl;
    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{
    for(int ict=0;ict<ct_Number;ict++){
      log( SIGNATURE, 
           "%d event(s), with a total of %d C.photons in CT %i (%s)\n", 
           ntshow, ntcph[ict],ict,GeometryName[ict] );
      datafile<<ntshow<<" event(s), with a total of "<<ntcph[ict]
              <<" C.photons in CT "<<ict<<" ("<<GeometryName[ict]<<")"<<endl;
      log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n", 
           ((float)ntrigger[ict]) / ((float)ntshow) * 100.0, ntrigger[ict], ntshow);
      datafile<<"Fraction of triggers: "<<((float)ntrigger[ict]) / ((float)ntshow) * 100.0<<" ("<<ntrigger[ict]<<" out of "<<ntshow<<" )"<<endl;
    }
  }

  // close files
  
  log( SIGNATURE, "Closing files\n" );
  
  if( ! Data_From_STDIN ){
    for(int ict=0;ict<ct_Number;ict++)
      fclose( inputfile[ict] );
  }
  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"
       << "O. Blanch, A. Moralejo, 2004\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_QE  
//                          
// @desc read QE data
//
// @date thu  5 17:59:57 CEST 2002
//------------------------------------------------------------
// @function

//!@{
void 
read_QE(char fname[256], int ict){
  ifstream qefile;
  char line[LINE_MAX_LENGTH];
  int i, j, icount;
  float qe;

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

  // try to open the file

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

  log("read_QE", "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[ict]);
      
      // allocate memory for the table of QEs
      
      QE[ict] = new float ** [ct_NPixels];

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

      for ( i=0; i<pointsQE[ict]; ++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 < ct_NPixels) && (i > -1) &&
         ((j-1) < pointsQE[ict])   && ((j-1) > -1) ) {
      QE[ict][i][0][j-1] = QElambda[j-1];
      QE[ict][i][1][j-1] = qe;
    }

    if ( i > ct_NPixels) break;

    icount++;

  }

  if(icount/pointsQE[ict] < ct_NPixels){
    error( "read_QE", "The quantum efficiency file is faulty\n  (found only %d pixels instead of %d).\n",
           icount/pointsQE[ict], ct_NPixels );
  }

  // close file

  qefile.close();

  // test QE

  for(icount=0; icount< ct_NPixels; icount++){
    for(i=0; i<pointsQE[ict]; i++){
      if( QE[ict][icount][0][i] < 100. || QE[ict][icount][0][i] > 1000. ||
          QE[ict][icount][1][i] < 0. || QE[ict][icount][1][i] > 100.){
        error( "read_QE", "The quantum efficiency file is faulty\n  pixel %d, point %d  is % f, %f\n",
               icount, i, QE[ict][icount][0][i], QE[ict][icount][1][i] );
      }
    }
  }

  // end

  log("read_QE", "Done.\n");
}
//!@}

//!-----------------------------------------------------------
// @name read_WC  
//                          
// @desc read WC data
//
// @date thu  5 17:59:57 CEST 2002
//------------------------------------------------------------
// @function

//!@{
void 
read_WC(void){
  ifstream wcfile;
  char line[LINE_MAX_LENGTH];
  int i;

  //------------------------------------------------------------
  // Read Light Collection data

  // try to open the file

  log("read_WC", "Opening the file \"%s\" . . .\n", WC_FILE);
  
  wcfile.open( WC_FILE );
  
  // if it is wrong or does not exist, exit
  
  if ( wcfile.bad() )
    error( "read_WC", "Cannot open \"%s\". Exiting.\n", WC_FILE );
  
  // read file

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

  // get line from the file

  do
    wcfile.getline(line, LINE_MAX_LENGTH);
  while (line[0] == '#');

  // get the number of datapoints 
    
  sscanf(line, "%d", &pointsWC);
    
  // allocate memory for the table of QEs
    
  WC = new float * [2];
  WC[0] = new float[pointsWC];
  WC[1] = new float[pointsWC];

  for ( i=0; i<pointsWC; ++i ) {
    wcfile.getline(line, LINE_MAX_LENGTH);
    sscanf(line, "%f %f", &WC[0][i], &WC[1][i]);
  }

  //
  // Now read info for outer pixels
  //
  WC_outer = new float * [2];
  WC_outer[0] = new float[pointsWC];
  WC_outer[1] = new float[pointsWC];

  do
    wcfile.getline(line, LINE_MAX_LENGTH);
  while (line[0] == '#');

  if (wcfile.eof())
    {
      log("read_WC", "ERROR. Missing data for outer pixels in file \"%s\"...\n",WC_FILE);
      log("read_WC", "EXITING camera\n");
      exit(-1);
    }

  sscanf(line, "%f %f", &WC_outer[0][0], &WC_outer[1][0]);
  for ( i=1; i<pointsWC; ++i ) {
    wcfile.getline(line, LINE_MAX_LENGTH);
    sscanf(line, "%f %f", &WC_outer[0][i], &WC_outer[1][i]);
  }

  // close file

  wcfile.close();

  // read

  log("read_WC", "Done.\n");
}
//!@}


//!-----------------------------------------------------------
// @name read_ascii  
//                          
// @desc read ascii configuration files used by the reflector
//
// @date tue dec 10 17:14:10 CET 2002
//------------------------------------------------------------
// @function

//!@{
void 
read_ascii(FILE *sp, MMcConfigRunHeader *config)
{
  Float_t radius = -1.0;
  Float_t focal = -1.0;
  Float_t point = -1.0;
  Float_t spot = -1.0;
  Float_t camwidth = -1.0;

  Int_t imir;
  Float_t f,sx,sy,x,y,z,thetan,phin,xn,yn,zn;
  Float_t dx,dy;
  Int_t nummir, numref;
  Float_t wav,ref;

  Char_t token[40];
  Char_t line[511];
  Char_t flag;

  while(1){
    if((flag=fgetc(sp))==EOF)
      break;
    if (flag == '\n')    // skip empty lines
      continue;

    fgets(&line[1],500,sp);
    line[0]=flag;

    if ( strstr(line, "# reflectivity file") == line ) {
      while (strstr(line, "# number of datapoints") != line)
        fgets(line,500,sp);

      fgets(line,500,sp);
      sscanf(line,"%i",&numref);

      TArrayF wavarray(numref);
      TArrayF refarray(numref);

      while (strstr(line, "# datapoints") != line)
        fgets(line,500,sp);

      for(int i=0; i<numref;i++){
        fgets(line,500,sp);
        if (line[0] == '#')
          {
            i--;
            continue;
          }
        sscanf(line,"%f %f",&wav,&ref);
        wavarray[i]=wav;
        refarray[i]=ref;
      }
      for (int j=0; j<nummir;j++){

        MGeomMirror &mirror = config->GetMirror(j);
        mirror.SetArraySize(numref);
        mirror.SetReflectivity(wavarray, refarray);
      }
      continue;
    }
    if (line[0]== '#')
      continue;
    if (strstr(line, "type") == line)
      continue;
    if (strstr(line, "focal_distance") == line){
      sscanf(line,"%s %f",token,&focal);
      continue;
    }
    if (strstr(line, "point_spread") == line){
      sscanf(line,"%s %f",token,&point);
      continue;
    }
    if (strstr(line, "black_spot") == line){
      sscanf(line,"%s %f",token,&spot);
      continue;
    }
    if (strstr(line, "camera_width") == line){
      sscanf(line,"%s %f",token,&camwidth);
      continue;
    }
    //
    // Skip obsolete magic.def entries:
    //
    if (strstr(line, "n_pixels") == line)
      continue;
    if (strstr(line, "pixel_width") == line)
      continue;
    if (strstr(line, "n_centralpixels") == line)
      continue;
    if (strstr(line, "n_gappixels") == line)
      continue;

    if (strstr(line, "n_mirrors") == line){
      sscanf(line,"%s %i",token,&nummir);
      config->InitSizeMirror(nummir);
      continue;
    }
    if (strstr(line, "r_mirror") == line){
      sscanf(line,"%s %f",token,&radius);
      continue;
    }
    if (strstr(line, "define_mirrors") == line){
      for(int i=0;i<nummir;i++){
        fgets(line,500,sp);
        sscanf(line,"%i %f %f %f %f %f %f %f %f %f %f %f",
               &imir,&f,&sx,&sy,&x,&y,&z,&thetan,&phin,&xn,&yn,&zn);
        config->AddMirror(i);
        MGeomMirror &mirror = config->GetMirror(i);
        mirror.SetMirrorContent(imir,f,sx,sy,x,y,z,thetan,phin,xn,yn,zn);
      }
      fgets(line,500,sp);

      while ( ! strstr(line, "axis deviation"))
        fgets(line,500,sp);

      for(int i=0;i<nummir;i++){
        fgets(line,500,sp);
        if (line[0] == '#')
          {
            i--;
            continue;
          }
        sscanf(line,"%f %f",&dx,&dy);
        MGeomMirror &mirror = config->GetMirror(i);
        mirror.SetMirrorDeviations(dx,dy);
      }      
      continue;
    }
  }
  config->SetMagicDef(radius, focal, point, spot, camwidth);
}


//!-----------------------------------------------------------
// @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;
  double ddyseg1, ddyseg2, ddyseg3, ddyseg4, ddyseg5;
  

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

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

int bpoint_is_in_pix(double dx, double dy, MGeomCam *pgeomcam)
{
    /* 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    numN  = pgeomcam->GetNumPixels();

    for (int i=0; i<numN; i++)
    {
      MGeomPix &pixel = (*pgeomcam)[i];
      const double b = pixel.GetD()/2;
      const double a = pixel.GetD()/sqrt3;
      
      const double xx = dx - pixel.GetX();
      const double yy = dy - pixel.GetY();

      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 (x+t*l, y+t*m, z+t*n) to point P(a,b,c)
// 
// We assume that vector (l, m, n) is normalized l^2+m^2+n^2 = 1
//------------------------------------------------------------

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

//------------------------------------------------------------
// @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 sign[20];                // auxiliary variable

  strcpy(sign, REFL_SIGNATURE_B);

  fread( (char *)sign, strlen(REFL_SIGNATURE_B), 1, infile);
  if (strcmp(sign, REFL_SIGNATURE_A) == 0){
    fread( (char *)sign, 1, 1, infile);
    return 4;
  }
  else if (strcmp(sign, REFL_SIGNATURE_B) == 0){
    fread( (char *)sign, 1, 1, infile);
    return 5;
  }
  else if (strcmp(sign, REFL_SIGNATURE_C) == 0){
    // An empty bin has been removed and therefore we do not need to rezd it.
    return 6;
  }
  else {
    cout << "ERROR: Signature of .rfl file is not correct\n";
    cout << '"' << sign << '"' << '\n';
    cout << "should be: " << REFL_SIGNATURE_A <<" or "<< REFL_SIGNATURE_B <<" or " << REFL_SIGNATURE_C <<" or "<< '\n';
    return(FALSE);
  }


}

//------------------------------------------------------------
// @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 size_rotated
//                                     
// @desc It rotates the NSB
//
// @date Tue Jan  7 17:09:25 CET 2003
// @function @code 
//------------------------------------------------------------

int size_rotated(
    float *rotated,
    float *nsb,
    float rho)
{
    int r=0; 
    float size_rotated;
    float Num_Pixels;
    float PixNum=0;
    float rho_pixel;
    int j=0;
    int k=0;

    for(int i=1; i<iMAXNUMPIX;i++){
// Substituir per Int_t radius[iMAXNUMPIX]={0,1,1,1,1,1,1,2,2,2,2, ...}
        Num_Pixels=0;
        for (int ii=1; ii<17;ii++){
            if (Num_Pixels >= i){
                r=ii-1;
                break;
            }
                
            
            
            if (ii<12){
                Num_Pixels+=ii*6;
                PixNum=6*ii;
                //      size_rotated = (360/PixNum);
                //rho_pixel = rho/size_rotated;  
            }   
            
            
            else
            {
                Num_Pixels+=(ii-6)*6;
                PixNum=(ii-6)*6;
                // size_rotated = (360/PixNum);
                // rho_pixel = rho/size_rotated;
            }
        }            
            
            
            
        size_rotated = (360/PixNum);
        rho_pixel = rho/size_rotated;
        
        //Buscar j i k que no canvin de r!!!
        
        j=i-int(rho_pixel);
        
        //cout<<"j inicial "<<j<<endl;
        int  MINPIXinner[13]= {0,1,7,19,37,61,91,127,169,217,271,331,397};
        int  MINPIXouter[5]={397,433,475,523,578};
            
        if( r<12 ){
            if (j < MINPIXinner[r])
            {
                j = i+6 * r - int(rho_pixel);
            }
            
            k=j-1;
            if (k < MINPIXinner[r])
            {
                k =  MINPIXinner[r+1]-1;
            }
        }
        if(  r > 11 && r < 16){    
            if (j < MINPIXouter[r-12])
            {
                j = i + (r - 6) * 6 - int(rho_pixel);
                
            }
            
            k=j-1;
            if ( k < MINPIXouter[r-12])
            {
                k = MINPIXouter[r-11] - 1;
            }
        }
        rotated[i]= (1 - ((rho_pixel)-floor(rho_pixel))) * nsb[j] + (rho_pixel-floor(rho_pixel)) * nsb[k];
    }
    return (int(rho_pixel));
    
}


//------------------------------------------------------------
// @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
                  class MGeomCam *camgeom, // 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, // number of photoelectrons in each pixel
                  int *incph,       // total number of cph within camera
                  float *tmin_ns,   // minimum arrival time of all phes
                  float *tmax_ns,    // maximum arrival time of all phes
                  int ict,      // Telescope that is being analised to get the right QE.
                  float mirror_fraction, // Fraction of total mirror really present
                  float fadc_jitter     // Random shift (max 1 slice) of pulses in 
                                        // FADC, to simulate FADC clock noise.
                   ){

  static uint i;
  static int k, ipixnum;
  static float cx, cy, wl, qe;
  static float cw;
  static MCCphoton photon;
  static float **qept;
  static char flag[SIZE_OF_FLAGS + 1];
  static float radius_mm;
  static UInt_t seed = (UInt_t)(get_seeds(0)*get_seeds(1));

  float time;
  float pmtjit;

  // reset variables

  for ( i=0; i<camgeom->GetNumPixels(); ++i )
    nphe[i] = 0.0;


  *itotnphe = 0;
  *incph = 0;

  radius_mm = camgeom->GetMaxRadius();

  TRandom random;
  random.SetSeed(seed);

  float C1, C2, C3, rho;

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

    photon.read(sp);
    photon.get_flag(flag);

    if (isA( flag, FLAG_END_OF_EVENT ))
      {
        fseek (sp, SIZE_OF_FLAGS-photon.mysize(), SEEK_CUR);
        break;  
      }

    //  Check if photon is inside trigger time range

    time = photon.get_t() ;

    if (time-*tmin_ns>TOTAL_TRIGGER_TIME) {
        continue;
    }

    //
    // Account for possibly missing mirrors, or lower reflectivity...
    // mirror_fraction is the fraction of the total mirror really 
    // working.
    //
    if (mirror_fraction < 1.)
        if (RandomNumber > mirror_fraction)
            continue;

    //
    // Pixelization
    //
        
    cx = photon.get_x();
    cy = photon.get_y(); 

    if (Spotsigma > 0.)
    {
        cx += random.Gaus(0.,Spot_x);
        cy += random.Gaus(0.,Spot_y);
    }
    if (missPointing > 0.)
    {
      // We should take intot acount the rotation of the FoV
      C1 = 0.48 * sin(Zenith) - 0.87 * cos(Zenith) * cos(Azimutal);
      C3 = (0.87 * cos(Zenith) - 0.48 * sin(Zenith) * cos(Azimutal));
      C2 = sqrt( sin(Zenith) * sin(Zenith) * sin(Azimutal) * sin(Azimutal) + C3 * C3 );
      rho = acos( C1/C2 );
      rho=(sin(Azimutal)<0 ? (2 * 3.14159 - rho) : rho);

      rho = rho*180/3.14159;

      cx += (missP_x*cos(rho)-missP_y*sin(rho))/(10*camgeom->GetConvMm2Deg());
      cy += (missP_x*sin(rho)+missP_y*cos(rho))/(10*camgeom->GetConvMm2Deg());
    }
    

    // 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)*10 > radius_mm ) ){ 
      continue;
            
    }

    ipixnum = bpoint_is_in_pix(cx*10, cy*10, camgeom);

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

    // AM changed meaning of incph: before it was all photons read from 
    // reflector file, now only those within a valid pixel:
    //
    // increase number of photons
    (*incph)++;

    //+++
    // QE simulation
    //---
    
    // set pointer to the QE table of the relevant pixel

    qept = (float **)QE[ict][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[ict]-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[ict]-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;

    //
    // Apply incident angular correction due to Light Guides, plexiglas,
    // 1st dynode collection efficiency, double crossings... etc. 
    // This information is contained in the file Data/LightCollection.dat
    //
    cw=photon.get_phi()*180./3.14159265;

    // find data point in the WC table (-> k)

    if ( camgeom->GetPixRatio(ipixnum) < 1.) // => Pixel is outer pixel
      {
        k = 0; 
        while (k < pointsWC-1 && WC_outer[0][k] < cw)
          k++;
        // correct the qe with WC data.
        qe = qe*lin_interpol(WC_outer[0][k-1], WC_outer[1][k-1], 
                             WC_outer[0][k],   WC_outer[1][k], cw);
      }

    else    // => Pixel is inner pixel
      {
        k = 0; 
        while (k < pointsWC-1 && WC[0][k] < cw)
          k++;
        // correct the qe with WC data.
        qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], cw);
      }


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

    //
    // PMT time jitter: gaussian, not negative (MTrigger::FillShow does
    // not accept negative times!)
    //
    do 
      pmtjit = random.Gaus(3.*pmt_jitter, pmt_jitter);
    while(pmtjit<0.);

    // store the new photoelectron

    fadc->Fill(ipixnum,
               (time-*tmin_ns) + pmtjit + fadc_jitter, 
               trigger->FillShow(ipixnum, time-*tmin_ns + pmtjit + fadc_jitter),
               !((*camgeom)[ipixnum].GetD()>(*camgeom)[0].GetD()));
    
    *itotnphe += 1;
  }

  seed = random.GetSeed(); // Get seed for next call
          
  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
                      MGeomCam *camgeom, // camera layout
                      float **rate_phepns,
                      // the product of this function:
                      // the NSB rates in phe/ns for each pixel 
                      int ict,
                      float mirror_fraction)
{
  uint i, j; // counters
  int k, ii; // counters

  static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
                                            WAVEBANDBOUND2,
                                            WAVEBANDBOUND3,
                                            WAVEBANDBOUND4,
                                            WAVEBANDBOUND5, 
                                            WAVEBANDBOUND6 };

  FILE *infile;    // the input file
  fpos_t fileposition; // marker on the input file
  static char cflag[SIZE_OF_FLAGS + 1]; // auxiliary variable
  static MCEventHeader evth; // the event header
  static MCEventHeader evth_2; // the event header
  static float nphe[iMAXNUMPIX];    // the number of photoelectrons in each pixel
  int reflector_file_version;
  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
  char flag_new[4]; 


  if (strlen(iname) == 0)
    {
      log( SIGNATURE, "No starfield input file has been provided.\n");
      return (0);
    }
  else // check if the starfield input file exists and open it
    {
      log(SIGNATURE, "Opening starfield input \"rfl\" file %s\n", iname );
      infile = fopen( iname, "r" );

      if ( infile == NULL )
        {
          log( SIGNATURE, "ERROR! Cannot open starfield input file: %s\n", iname );
          exit(-1);
        }
    }


  // get signature, and check it
  
  if((reflector_file_version=check_reflector_file( infile ))==FALSE){
    exit(1);
  }

  // Instance of MTrigger and MFadc; needed here only as dummies for 
  // a call to produce_phes (see below).
  // 15/09/2004, A. Moralejo, changed "trigger" and "flashadc" to 
  // pointers, the former static allocation caused memory problems and
  // segmentation fault in some systems.

  MTrigger* trigger = new MTrigger(camgeom->GetNumPixels(),camgeom,
                                   Trigger_gate_length, 
                                   Trigger_overlaping_time, 
                                   Trigger_response_ampl, 
                                   Trigger_response_fwhm);

  MFadc* flashadc = new MFadc(camgeom->GetNumPixels(),
                              FADC_shape,
                              FADC_response_integ,FADC_response_fwhm,
                              FADC_shape_out,
                              FADC_resp_integ_out,FADC_resp_fwhm_out,
                              get_trig_delay(),
                              FADC_slices_per_ns,
                              FADC_slices_written);

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

      if (reflector_file_version > 5){

        fread( flag_new, 4, 1, infile );

        if(!isA( flag_new, FLAG_START_OF_HEADER)){
          
          //  We break the main loop
          cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
          cout<<"         We assume no Star Light"<<endl;
          return(0);
        }
        
        Float_t flag_temp[SIZE_OF_MCRUNHEADER];
        
        fread( flag_temp, (SIZE_OF_MCRUNHEADER)*sizeof(float), 1, infile );
        
      }
      
      fread( cflag, SIZE_OF_FLAGS, 1, infile );
      
      if( isA( cflag, FLAG_START_OF_EVENT   )){ // there is a event
        fread( flag_new, 4, 1, infile );

        if(!isA( flag_new, FLAG_EVENT_HEADER)){

          //  We break while events loop
          cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<endl;
          cout<<"         We assume no light from Stars"<<endl;
          return(0);      
        }

        // get MCEventHeader
        
        if (reflector_file_version<6)
          //      fread( (char*)&evth, evth.mysize(), 1, infile );
          evth.read(infile);
        else
          //      fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
          evth_2.read(infile);

        if (reflector_file_version<6)
          integtime_ns = evth.get_energy();
        else
          integtime_ns = evth_2.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<camgeom->GetNumPixels(); 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 (photons from starfield)

            k = produce_phes( infile,
                              camgeom,
                              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,
                              0,
                              mirror_fraction,
                              0.); 


            if( k != 0 ){ // non-zero returnvalue means error
              cout << "Exiting.\n";
              exit(1);
            }
          
            for(j = 0; j<camgeom->GetNumPixels(); 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);
}



//-------------------------------------------------------------
//
// Function DoCalibration. A. Moralejo October 2004
//
// Generates calibration events similar to those in the real 
// calibration runs.
//
//-------------------------------------------------------------

int DoCalibration(MFadc **Fadc_CT, MTrigger **Trigger_CT,
                  TObjArray camgeom,
                  float *nsb_trigresp, float *nsb_fadcresp,
                  MLons *lons, MLons *lons_outer,
                  float **nsb_phepns, int addElecNoise,
                  TTree *EvtTree, MRawEvtHeader **EvtHeader,
                  MMcEvt **McEvt, MRawEvtData** EvtData)
{

  int nevent = 0;
  int maxevents;
  int lons_return;

  int *ntotphe = new int[ct_Number];
  int *ntotphot = new int[ct_Number];
  float *nphe_from_nsb = new float[ct_Number];
  float *pulsetime = new float[ct_Number];
  
  float lambda, sigma_lambda, phot_per_pix, sigma_time;
  UInt_t *first_pixel, *last_pixel;
  int selected_pixel;

  TRandom random;
  random.SetSeed((UInt_t)(get_seeds(0)*get_seeds(0)));

  for (int ict = 0; ict < ct_Number; ict++)
    pulsetime[ict] = 0.;

  // Get parameters of the calibration:

  get_calibration_properties( &lambda, &sigma_lambda, &phot_per_pix, 
                              &sigma_time, &maxevents, &selected_pixel);

  first_pixel = new UInt_t[ct_Number];
  last_pixel  = new UInt_t[ct_Number];

  for (int ict = 0; ict < ct_Number; ict++)
    {
      if (selected_pixel < 0)
        {
          first_pixel[ict] = 0;
          last_pixel[ict]  = ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels()-1;
        }
      else
        first_pixel[ict] = last_pixel[ict] = selected_pixel;
    }


  TArrayC *fadcValues;     //@< the analog Fadc High gain signal for pixels
  TArrayC *fadcValuesLow;  //@< the analog Fadc Low gain signal for pixels

  // allocate space for FADC info
  fadcValues    =  new TArrayC(FADC_slices_written);
  fadcValuesLow =  new TArrayC(FADC_slices_written);


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

  for (int calevent = 0; calevent < maxevents; calevent++) 
    { 
      //
      // Clear Trigger and Fadc
      //
      for(int ict = 0; ict < ct_Number; ict++)
        {
          Trigger_CT[ict]->Reset() ; 
          Trigger_CT[ict]->ClearFirst();
          Trigger_CT[ict]->ClearZero();
          Fadc_CT[ict]->Reset() ; 

          ntotphe[ict] = ntotphot[ict] = 0;
          nphe_from_nsb[ict] = 0.;

        }

      nevent++;

      if((nevent+1)%100 == 1)
        log(SIGNATURE, "Event %d\n", nevent);
        
      // Produce the photoelectrons
      for(int ict = 0; ict < ct_Number; ict++)
        {
          // Obtain the FADC jitter of 1 FADC slice. This is a time to be added to the 
          // time of all photons in an event, before digitalization of the signal. It is
          // therefore the same time shift for all pixels in a CT. 

          float fadc_jitter = 
            (1./Fadc_CT[ict]->GetFadcSlicesPerNanosec()) * random.Uniform(0., 1.); //ns


          produce_calib_phes( (MGeomCam*)(camgeom.UncheckedAt(ict)),
                              Trigger_CT[ict],
                              Fadc_CT[ict],
                              &(ntotphe[ict]),
                              pheinpix,
                              &(ntotphot[ict]),
                              ict,
                              lambda,
                              sigma_lambda,
                              phot_per_pix,
                              sigma_time,
                              selected_pixel,
                              fadc_jitter
                              );

          pulsetime[ict] = fadc_jitter; // Keep value for writing it to output.
        }

      // NSB simulation

      if(lons && lons_outer)
        {
          //  Fill trigger and fadc response in the trigger class from the NSB database
          
          for (int ict = 0; ict < ct_Number; ict++)
            {
              for (UInt_t ui = first_pixel[ict]; ui <= last_pixel[ict]; ui++)
                {
                  nphe_from_nsb[ict] += nsb_phepns[ict][ui];

                  if (nsb_phepns[ict][ui] > 0.0)
                    {
                      if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD() > 
                         (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD())
                        lons_return = lons_outer->GetResponse(nsb_phepns[ict][ui],0.01,
                                                              & nsb_trigresp[0],
                                                              & nsb_fadcresp[0]);
                      else
                        lons_return = lons->GetResponse(nsb_phepns[ict][ui],0.01,
                                                        & nsb_trigresp[0],
                                                        & nsb_fadcresp[0]);

                      if (lons_return == 0)
                        {
                          cout << "Exiting.\n";
                          exit(1);
                        }

                      Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
                    }
                }
            }
          
        }// end if(simulateNSB) ...


      //++++++++++++++++++++++++++++++++++++++++++++++++++
      // at this point we have a camera full of
      // ph.e.s
      //--------------------------------------------------
          
      //   now the noise of the electronic 
      //   (preamps, optical transmission,..)  is introduced. 
      //   
      
      for (int ict = 0; ict < ct_Number; ict++) 
        {
          if (addElecNoise)
            Fadc_CT[ict]->ElecNoise();

          //  now a shift in the fadc signal due to the pedestals is 
          //  introduced
          //  This is done inside the class MFadc by the method Pedestals

          Fadc_CT[ict]->Pedestals();

          // Set the trigger time. The 3 ns are such that the calibration pulses
          // appear roughly at the same position as for the case of real data.
          // If you want to shift the pulse position, do not change this value here,
          // use the option trigger_delay in the input card instead.
          // The additional value 3*sigma_time makes that the pulse maximum is, in 
          // average, in the same position no matter of the time width of the pulse
          // (see that, in produce_calib_phes(...), in order to avoid negative times 
          // we shift the time of the photons by this same amount!)

          Fadc_CT[ict]->TriggeredFadc(3.+3*sigma_time);

          // Add the "digital noise": electronic noise intrinsic to the FADC and which
          // therefore is not scaled down in the low gain slices!

          Fadc_CT[ict]->DigitalNoise();

          //
          // Fill the event header information
          //
          EvtHeader[ict]->FillHeader((UInt_t) calevent, 0);

          // Set flag to indicate this is a calibration or a pedestal run:
          if (phot_per_pix > 1.e-3)
            {
              EvtHeader[ict]->SetTriggerPattern((UInt_t)(MTriggerPattern::kCalibration |
              MTriggerPattern::kTriggerLvl1));
              //
              // FIXME! For now color and intensity of the pulser is fixed!
              EvtHeader[ict]->SetCalibrationPattern((UInt_t)((BIT(11) | BIT(12))<<16));
            }
          else  // 0 cal. photons per pixel ==> pedestal run
            EvtHeader[ict]->SetTriggerPattern((UInt_t)MTriggerPattern::kPedestal);


          McEvt[ict]->Fill( 0, 0, 0., -1.0, -1.0, -1.0, 0., 0., 0., 0., 0.,
                            0., 0., 0., 0., 0., 0., 0., 0.,
                            0., 0., 0., 0., 0, 0, 0,
                            (UInt_t) ntotphot[ict], 
                            (UInt_t) ntotphe[ict],
                            (UInt_t) nphe_from_nsb[ict]+ ntotphe[ict],
                            0., 0., 0., pulsetime[ict]);

          //
          // Fill the FADC information
          //
          for(UInt_t ui = first_pixel[ict]; ui <= last_pixel[ict]; ui++)
            {
              for (int jslice = 0; jslice < FADC_slices_written; jslice++)
                {
                  fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(ui, jslice),jslice);
                  fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(ui,jslice),jslice);
                }
              EvtData[ict]->AddPixel(ui,fadcValues,0);
              EvtData[ict]->AddPixel(ui,fadcValuesLow,kTRUE);
            }
        }

      EvtTree->Fill(); 

      //    clear all
      for(int ict=0;ict<ct_Number;ict++)
        {
          EvtHeader[ict]->Clear() ;
          EvtData[ict]->ResetPixels(0,0);
          McEvt[ict]->Clear() ; 
        }
    }

  return(0);
}


//------------------------------------------------------------
// 
// Function produce_calib_phes, A. Moralejo Oct 2004
//    
//------------------------------------------------------------

int produce_calib_phes( MGeomCam *camgeom,       // The camera layout
                        MTrigger *trigger,
                        MFadc    *fadc,
                        int      *itotnphe,      // total number of produced photoelectrons             
                        float    *nphe,          // number of photoelectrons in each pixel
                        int      *nphot,         // total number of photons in all pixels
                        int       ict,           // Telescope that is being analised to get the right QE.
                        float     lambda,        // Mean wavelength of light in nm
                        float     sigma_lambda,  // Sigma of wavelengtgaussian spread
                        float     phot_per_pix,  // Average # of photons per inner pixel
                        float     sigma_time,    // Sigma of time spread of photons
                        int       selected_pixel,// if >= 0, only this pixel is used! 
                        float     fadc_jitter    // Random shift (max 1 slice) of pulses in 
                                                 // FADC, to simulate FADC clock noise.
                        )
{

  int   ipixnum;
  float cx, cy, wl, time, phi, phi_deg, qe;
  float **qept;
  float radius_mm, focal_dist_mm;
  int   total_photons;
  float pmtjit;

  static UInt_t seed = (UInt_t)(get_seeds(0)*get_seeds(1));

  // reset variables

  for ( uint i = 0; i < camgeom->GetNumPixels(); i++ )
    nphe[i] = 0.0;

  *itotnphe = 0;
  *nphot = 0;

  TRandom random;
  random.SetSeed(seed);

  //
  // Create photons and "map" them into the pixels
  //
  
  // Maximum radius of camera:
  radius_mm = camgeom->GetMaxRadius();

  // Distance from center of mirror dish to camera plane:
  focal_dist_mm = camgeom->GetCameraDist()*1000;

  // Cosine of the angle between telescope axis and line from center of mirror
  // dish to the edge of the camera:
  float cos_epsilon_max = cos(atan2(radius_mm, focal_dist_mm));


  // Calculate total number of photons to be produced.
  if (selected_pixel < 0)
    total_photons = (int) (phot_per_pix * 3.14159265 * radius_mm * radius_mm / 
                           (*camgeom)[0].GetA());
  else
    total_photons = (int) (phot_per_pix * 
                           (*camgeom)[selected_pixel].GetA() / (*camgeom)[0].GetA());

  // loop over the photons
  
  for (int iph = 0; iph < total_photons; iph++) 
    {
      //
      // Simulate arrival times of photons to camera plane. We do not simulate the small
      // delays between pixels due to the different distances to the pulser.
      //
      // We do not want negative times, so center the gaussian at 3 sigma 
      // and reject negative values:

      do
        time = random.Gaus(3*sigma_time, sigma_time);
      while (time < 0.);

      // wavelength
      wl = random.Gaus(lambda, sigma_lambda);

      if( (wl > WAVEBANDBOUND6) || (wl < WAVEBANDBOUND1))
        continue;

      if (selected_pixel < 0)
        {
          // Obtain photon coordinates on the camera. We assume a point source of light placed
          // in the center of the mirror dish.

          // polar angle 
          float psi = RandomNumber * 2 * 3.14159265;
          // angle between the telescope axis and the photon trajectory.
          float epsilon = acos(1.-RandomNumber*(1.-cos_epsilon_max));
          float tanepsilon = tan(epsilon);

          cx = focal_dist_mm*tanepsilon*cos(psi); // mm
          cy = focal_dist_mm*tanepsilon*sin(psi); // mm
 
          if (sqrt(cx*cx + cy*cy) > radius_mm )
            continue;

          // Angle between photon trajectory and camera plane:
          phi = 3.14159265/2.-epsilon; // rad

          //
          // Pixelization
          //
          ipixnum = bpoint_is_in_pix(cx, cy, camgeom);

          // -1 = the photon is in none of the pixels
          //  0 = the phton is in the central pixel, which is not used
          if (ipixnum==-1 || ipixnum==0)
            continue;
        }
      else
        {
          // Angle between photon trajectory and camera plane:
          phi = atan2( focal_dist_mm, 
                       sqrt( (*camgeom)[selected_pixel].GetX()*(*camgeom)[selected_pixel].GetX()+
                             (*camgeom)[selected_pixel].GetY()*(*camgeom)[selected_pixel].GetY()));

          ipixnum = selected_pixel;
        }


      // increase number of photons within pixels
      *nphot += 1;

      //
      // QE simulation
      //    
      // set pointer to the QE table of the relevant pixel
      
      qept = (float **)QE[ict][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[ict]-1]))
        continue;
            

      // find data point in the QE table (-> k)
      int k = 1; // start at 1 because the condition was already tested for 0
      while (k < pointsQE[ict]-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;


      //
      // Apply incident angular correction due to Light Guides, plexiglas,
      // 1st dynode collection efficiency, double crossings... etc. 
      // This information is contained in the file Data/LightCollection.dat,
      // and has been read into the array WC (which stands for "Winston Cones")
      //
      phi_deg = phi*180./3.14159265;

      // find data point in the WC table (-> k)

      if ( camgeom->GetPixRatio(ipixnum) < 1.) // => Pixel is outer pixel
        {
          k = 0; 
          while (k < pointsWC-1 && WC_outer[0][k] < phi_deg)
            k++;
          // correct the qe with WC data.
          qe = qe*lin_interpol(WC_outer[0][k-1], WC_outer[1][k-1], 
                               WC_outer[0][k],   WC_outer[1][k], phi_deg);
        }

      else    // => Pixel is inner pixel
        {
          k = 0; 
          while (k < pointsWC-1 && WC[0][k] < phi_deg)
            k++;
          // correct the qe with WC data.
          qe = qe*lin_interpol(WC[0][k-1], WC[1][k-1], WC[0][k], WC[1][k], phi_deg);
        }


      // if random > quantum efficiency, reject it

      if ( (RandomNumber) > qe ) {
        continue;
      }
          
      //
      // The photon has produced a photo electron
      //

      // increment the number of photoelectrons in the relevant pixel
          
      nphe[ipixnum] += 1.;

      //
      // PMT time jitter: gaussian, not negative (MTrigger::FillShow does
      // not accept negative times!)
      //
      do 
        pmtjit = random.Gaus(3.*pmt_jitter, pmt_jitter);
      while(pmtjit<0.);

      // store the new photoelectron

      fadc->Fill(ipixnum, time + pmtjit + fadc_jitter, 
                 trigger->FillShow(ipixnum, time + pmtjit + fadc_jitter),
                 !((*camgeom)[ipixnum].GetD()>(*camgeom)[0].GetD()));
    
      *itotnphe += 1;
    }

  seed = random.GetSeed(); // Get seed for next call

  return(0);

}

// @endcode


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

//!@{
//
// $Log: not supported by cvs2svn $
// Revision 1.88  2005/02/11 20:00:01  moralejo
//
// Added to output container "MMcEvtBasic" with the most important MC
// parameters, to be kept for all events (triggered or not) through
// the whole analysis chain in order to allow the calculation of effective
// areas.
//
// Updated version of MRawRunHeader from 4 to 5. This means the trigger and
// calibration patterns are correctly set and can be decoded.
//
// Revision 1.87  2005/02/10 19:28:10  moralejo
//
// Substituted input card option "write_all_events" by
// "write_all_event_headers". If set, this will make camera to write out
// the event headers of all the processed events, no matter whether they
// have triggered or not.
//
// Revision 1.86  2005/02/10 14:49:21  moralejo
//
// Added setting of trigger flags in MRawEvtHeader::fTriggerPattern:
// Lvl1 flag for normal shower runs, Calibration flag for calibration
// runs and Pedestal flag for "calibration runs" in which the number
// of photons per pixel is set to 0.
//
// Revision 1.85  2005/02/10 12:00:32  moralejo
//
//  Changed call to EvtHeader[ict]->FillHeader in line 2535: second argument
//  was 20, and I set it now to 0. No idea why that was 20: like that, it was
//  setting a CL trigger pattern on all events (???!!!)
//
// Revision 1.84  2004/12/15 01:56:39  moralejo
//
//  Added input card option pmt_jitter_ns to simulate the time jitter of
//  PMTs. The input parameter is the sigma of a gaussian, which by
//  default is sigma=0.25 ns. This jitter is applied to each phe-
//  independently. We have not applied this to the NSB noise, since the
//  arrival time of NSB photons is random and nothing would change.
//
// Revision 1.83  2004/11/17 11:43:13  moralejo
//
//  Made the necessary changes in starresponse to account for the new option
//  to switch off gain fluctuations for the NSB noise database generation.
//  The option in the Star Response card is "gain_fluctuations_off". The
//  version number of the NSB database has been updated to 1004 (see
//  MStarLight.hxx), now including information on whether the gain
//  fluctuations were on or off when the NSB database was generated.
//
// Revision 1.82  2004/11/17 11:34:49  moralejo
//
//  Added input card command  "noise_gain_fluctuations_off", to disable the
//  PMT gain fluctuations also for the NSB noise (just for tests). Added
//  a flag to MMcFadcHeader and MMcTrigHeader about this. Also copied the flag
//  fGainFluctuations of MMcFadcHeader to MMcTrigHeader, since the gain
//  fluctuations affect both the trigger and the signal in the FADC.
//
// Revision 1.81  2004/11/04 22:00:51  moralejo
//
//  Removed unused variables fTelesTheta and fTelesPhi from MMcRunHeader. They
//  were not useful because telescope orientation may change from event to
//  event. Added fMirrorFraction to the MMcConfigRunHeader container in the
//  camera output.
//
// Revision 1.80  2004/10/26 19:21:05  moralejo
//
//  Added fFadcTimeJitter to root output, in container MMcEvt. Added 
//  also fGainFluctuations boolean flag to MMcFadcHeader, to keep track of
//  whether PMT gain fluctuations are simulated or not.
//
// Revision 1.79  2004/10/26 14:02:32  moralejo
//
// Added option to switch off gain fluctuations (gain_fluctuations_off)
//
// Revision 1.78  2004/10/21 17:44:07  moralejo
//
//  Fixed error recently introduced in MLons::GetResponse. The NSB database
//  was not correclt ycopied to the FADC when the randomly selected time was
//  too close to the end of the database. This happened about 2% of times and
//  produced some repetitive noise peaks in the FADC.
//  Changed in camera.cxx the arguments of lons.SetSeed y lons_outer.SetSeed
//
// Revision 1.77  2004/10/19 10:35:05  moralejo
// *** empty log message ***
//
// Revision 1.76  2004/10/14 16:56:43  moralejo
//
// - Added calibration_run option to produce calibration MC files.
//
// - Added jitter of pulse position +- 0.5 slices due to FADC clock noise.
//
// Revision 1.75  2004/10/14 12:55:02  moralejo
// *** empty log message ***
//
// Revision 1.74  2004/10/13 17:05:05  moralejo
// *** empty log message ***
//
// Revision 1.72  2004/10/12 13:39:34  moralejo
//
// Lots of changes intended to make it possible to select the FADC sampling
// frequency from the input card, through the command fadc_GHz. The most
// important ones are the following:
//
//  - Removed FADC_SLICES de Mars/mmc/Mdefine.h
//    Already defined in MFadcDefine.h!
//
//  - Replaced fixed numbers in array dimensions of starresponse.cxx
//
//  - Added in MFadc.cxx and MFadc.hxx (Float_t) casts to initializations
//    of single phe response array
//
//  - IMPORTANT: Fixed MFadc::GetResponse  -> The returned single phe response
//    had only RESPONSE_SLICES (which is actually for the trigger branch),
//    whereas it should have RESPONSE_SLICES_MFADC. Fixed the same confusion
//    in other two points of the code (filling of the FADC for the signal),
//    in Fill() and FillOuter().
//
//  - RESPONSE_SLICES_MFADC is now eliminated, since this quantity is now
//    decided by MFadc depending on other parameters.
//
//  - Fixed problem in starresponse.cxx due to which the histograms fadcresp
//    and fadcbase in the root output were actually identical.
//
//  - Changed starresponse.cxx such that above 1 phe/ns/pixel the precision
//    of the database is forced to be 0.1, and below 1, to 0.01 phe/ns/pixel
//
//  - In Camera/creadparam.cxx: Now unkown tokens cause the program to stop,
//    instead of being simply skipped as it was until now.
//
//  - Fixed error in MStarLight::StoreHisto. TH1::SetBinContent uses bin numbers
//    from 1 to nbins (in old code 0 to nbins-1 was assumed).
//
//  - In MLons.cxx: use memcpy to copy pieces of the database into the FADC and
//    trigger branches, to (hopefully) speed up execution. For this I had to add
//    2 getter functions in MStarLight.hxx
//
//  - Everywhere: changed the shape parameter for trigger and FADC response to
//    be an Int. Changed version of NSB database from 1002 to 1003.
//
//  - In MTrigger.cxx, changed all initializations of SlicesFirst and
//    SlicesSecond to 0 (instead of -50 as it was before). This controls the
//    time of trigger. If no trigger happened (like when making pedestal files)
//    the time was negative and the array index used to retrieve the noise from
//    the database was negative, resulting in "discontinuities" in the noise
//    (half-photoelectrons for instance).
//
//  - In MStarLight changed fBinsTrig  and fBinsFadc from Float_t to Int_t
//
//  - Replace WIDTH_FADC_TIMESLICE by FADC_SLICES_PER_NSEC (which is the
//    inverse of the former)
//
//  - Replaced SLICES_PER_NSEC by TRIG_SLICES_PER_NSEC
//
//  - TRIGBINS eliminated. It depends on other two values
//        TRIGBINS = TIMERANGE*TRIG_SLICES_PER_NSEC
//
//  - FADCBINS eliminated. It depends on other two values
//        FADCBINS = TIMERANGE*FADC_SLICES_PER_NSEC
//
//  - MTriggerDefine.h  Changed RESPONSE_SLICES to RESPONSE_SLICES_TRIG
//
//  - Added to the MLons constructor an argument regarding the FADC sampling
//    frequency
//
//  - MFadc: now the number of response slices for the FADC simulation is
//    decided by the program according to the other parameters.
//
//  - MStarLight: Added arguments (fadc_slices_per_ns, response_slices_fadc)
//    to constructor.
//
//  - creadparam.cxx, camera.cxx   Changed default value of digital_noise to 0.
//
// Revision 1.71  2004/09/17 09:20:52  moralejo
//
// Updated some calls to current version of Mars:
//
// - EvtData[i]->InitRead(RunHeader)  instead of  EvtData[i]->Init(RunHeader);
// - MRawRunHeader::kMagicNumber instead of just kMagicNumber
// - EvtData[i]->ResetPixels (0, 0) instead of EvtData[i]->DeletePixels();
//
// Revision 1.70  2004/09/16 15:23:12  moralejo
//
// Changed "flashadc" and "trigger" in procedure produce_nsbrates from
// objects to pointers (followed by dynamical allocation). This is only
// to avoid memory problems (-> segmentation fault) in some systems.
// Introduced missing initialization to 0 of *itotnphe in produce_phes.
// Now the number of phes produced by stars shown on the screen make sense.
//
//
// Revision 1.69  2004/03/30
// Changed calculation of MMcFadcHeader.fPedesSigmaHigh and 
// MMcFadcHeader.fPedesSigmaLow to do as in real data (see comments in
// code). Changed meaning of MMcFadcHeader.fAmplFadc and fAmplFadcOuter,
// from amplitude to integral of single photoelectron pulse in FADC 
// counts. Added possibility to choose a realistic pulse shaped (as 
// measured using pulpo). Changed file Data/lightguides.dat by 
// Data/LightCollection.dat, intended to contain the information on 
// light collection efficiency regarding Winston cones, plexiglas, double 
// PMT crossing and colection efficiency of 1st dynode of PMT. Now the
// information for inner and outer pixels can be different, since in the
// LightCollection.dat file they are set independently.
//
// Revision 1.68  2004/02/06 17:39:24  blanch
// Compiling with root 3.05 and updating MARS files.
//
// Revision 1.67  2004/01/30 09:51:18  blanch
// [Changes mainly done by A. Moralejo]
//
// Several minnor changes have been done. For instance, some name of the
// variables have been modified to a more self-explained name and
// modifications while reading the asciis files at the end of the reflector file
// has been introduced.
//
// The possibilty to enlarge the point spread function has been introduced
// in order to be able to simualte the current data.
//
// All pixels are always written.
//
// Revision 1.65  2003/10/26 19:43:00  blanch
// - The screen output information has been improved to prevent some
// non-desired running conditions just looking at the output screen.
// - One MMcEvt for each Telsecope is stored in the output file.
// - 500 empty events are simualted to get a more precise estimation of the
// pedestal Sigma for each pixel.
//
// Revision 1.64  2003/10/21 07:42:50  blanch
// A factor 2.35 to transform the fwhm into the sigma of gaussian was missing
// in the storing of FADC single hpe pulse determination.
//
// Revision 1.63  2003/10/17 19:38:31  blanch
// Now the camera program will stop if a undefined Geometry is required.
// The NSB is internally scaled for any camera geometry and qe.
// The seeds in the input card are used to initilize the random numbers.
// The Amplitud stored in the MMcFadcHeader is the amplitud of the sphe reponse.
// The Pedestal rms is simulated in an artificial empty event.
//
// Revision 1.62  2003/09/26 11:25:07  blanch
// Modification to be able to read MGeomCam branch for any Geometry.
//
// Revision 1.61  2003/09/25 17:09:20  blanch
// Bug on the number of phe from diffuse NSB fixed.
//
// Revision 1.60  2003/09/23 16:50:55  blanch
// WE do not read ct_file anymore since all Telescope information is
// in the reflector or in MGeomCam.
//
// Revision 1.58  2003/09/15 09:59:53  blanch
// The concept of the camera prgoram has not changed but this version has
// quite a lot of changes to allow several Camera geometries as well as
// multitelescope.
//
// It is suposed to be the first working code for camera 0.7.
//
// Revision 1.57  2003/07/17 18:02:46  blanch
// Several new features introduced as well as fixed bugs
//
//      - 1/100 events printed out
//      - Low gain implemented
//      - Different response for outer and inner pixels
//      - Some warnings removed
//      - pedestal and qe file from inpuit card
//      - Faster electronic simulation
//
// Revision 1.55  2003/02/12 12:22:10  blanch
// *** empty log message ***
//
// Revision 1.54  2003/02/12 11:55:01  blanch
// *** empty log message ***
//
// Revision 1.53  2003/01/23 18:35:21  blanch
// *** empty log message ***
//
// Revision 1.52  2003/01/20 17:19:20  blanch
// It fills the MMcCorsikaRun.
//
// Revision 1.51  2003/01/14 13:40:17  blanch
// MRawRunHeader::fNumEvents has been filled with the number of events in
// this file.
// Problems in fImpact computation have been solved.
// Option to set a dc value to rise the discriminator threshold has been added.
//
// Revision 1.50  2003/01/07 16:33:31  blanch
// Star Field Rotation has been implemented by Raul Orduna. Now there is a
// rotation for each shower. It is done by a non enter pixel rotation assuming
// a circular symetry of the camera. It is not exact but it is accurate enough and
// much faster.
//
// Revision 1.49  2002/12/13 10:04:07  blanch
// *** empty log message ***
//
// Revision 1.48  2002/12/12 17:40:50  blanch
// *** empty log message ***
//
// Revision 1.47  2002/12/10 17:19:31  blanch
// *** empty log message ***
//
// Revision 1.46  2002/11/08 17:51:00  blanch
// *** empty log message ***
//
// Revision 1.45  2002/10/29 17:15:27  blanch
// Reading several reflector versions.
//
// Revision 1.44  2002/10/18 16:53:03  blanch
// Modification to read several reflector files.
//
// Revision 1.43  2002/09/13 10:53:39  blanch
// Minor change to remove some undisired comments.
//
// Revision 1.42  2002/09/09 16:00:49  blanch
// Statement has been included to avoid writting to disk MParContainer and MArray.
// It has also been added the effect of the WC, the actual values must be added,
// once they are measured.
//
// Revision 1.41  2002/09/04 09:57:42  blanch
// Modifications done to use MGeomCam from MARS.
//
// Revision 1.40  2002/07/16 16:15:22  blanch
// A first implementation of the Star field rotation has been done.
//
// Revision 1.39  2002/04/27 10:48:39  blanch
// Some unused varibles have been removed.
//
// Revision 1.38  2002/03/18 18:44:29  blanch
// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
//
// Revision 1.37  2002/03/18 16:42:20  blanch
// The data member fProductionSite of the MMcRunHeader has been set to 0,
// which means that the prodution site is unknown.
//
// Revision 1.35  2002/03/15 15:15:52  blanch
// Several modifications to simulate the actual trigger zone.
//
// Revision 1.34  2002/03/13 18:13:56  blanch
// Some changes to fill correctly the new format of MMcRunHeader.
//
// Revision 1.33  2002/03/04 17:21:48  blanch
// Small and not important changes.
//
// Revision 1.32  2002/02/28 15:04:52  blanch
// A small back has been solved. Before, while not using the option
// writte_all_images, not all triggered showers were stored. Now it is solved.
// For that it is important that the less restrictive trigger option is
// checked first.
// A new facility has been introduced and now one can choose the step size in
// trigger loop mode for the discriminator threshold.
// The close-compact topology for less than 3 pixels does not make sense. Before
// the program was ignoring that, now it switch to simple neighbour condition.
//
// 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