Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 3044)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 5072)
@@ -21,7 +21,7 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.68 $
-// $Author: blanch $ 
-// $Date: 2004-02-06 17:39:24 $
+// $Revision: 1.69 $
+// $Author: moralejo $ 
+// $Date: 2004-09-16 15:16:34 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -84,4 +84,5 @@
 #include "MGeomCamCT1Daniel.h"
 #include "MGeomPix.h"
+#include "MGeomCorsikaCT.h"
 
 /*!@" 
@@ -151,7 +152,4 @@
 static int ct_Number;        
 
-//@: Camera geometries
-static int ct_Geometry;        
-
 //@: list of showers to be skipped
 static int *Skip;
@@ -200,9 +198,12 @@
 
 //@: Properties of the FADC
-static float FADC_response_ampl = MFADC_RESPONSE_INTEGRAL;
+static int FADC_shape = 0;
+static float FADC_response_integ = MFADC_RESPONSE_INTEGRAL;
 static float FADC_response_fwhm = MFADC_RESPONSE_FWHM;
-static float FADC_resp_ampl_out = MFADC_RESPONSE_INTEGRAL;
+static int 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_noise = 2.0;
+static float FADC_noise_inner = 2.0;
+static float FADC_noise_outer = 2.0;
 static float DIGITAL_noise = 2.0;
 static float FADC_high2low = HIGH2LOWGAIN;
@@ -210,6 +211,6 @@
 //@: Trigger conditions for a single trigger mode
 static float **qThreshold;  
-static int Trigger_multiplicity[100];
-static int Trigger_topology[100];
+static int Trigger_multiplicity[MAX_NUMBER_OF_CTS];
+static int Trigger_topology[MAX_NUMBER_OF_CTS];
 
 //@: Upper and lower edges of the trigger loop
@@ -222,4 +223,19 @@
 static int Trigger_loop_utop = 2;
 
+//@: Direction of each shower
+static  float Zenith = 0.0;
+static  float Azimutal = 90.0;
+
+//@: Miss Pointing 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;
+
+  
 //!@}
 
@@ -259,11 +275,12 @@
 
 //@: number of datapoints for the QE curve
-static int pointsQE[100];
+static int pointsQE[MAX_NUMBER_OF_CTS];
 
 //@: table of QE
 static float *QElambda;
 
-//@: table of lightguide = WC;
+//@: table of lightguide = WC;  WC_outer for outer pixels.
 static float **WC;
+static float **WC_outer;
 
 //@: number of datapoints for the WC curve
@@ -313,5 +330,5 @@
   //@'
 
-  char inname_CT[100][256];   //@< array of names for each CT input file
+  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
@@ -328,40 +345,46 @@
 
   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[100][50];            //@< Name of MGeomCam for each CT
-  int GeometryCamera[100];           //@< Identification of MGeomCam for each CT
-  int TriggerPixels[100];           //@< Number of pixels in the trigger region for each CT
+  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[100];            //@< stream for the input 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[100];       //@< Event Header class (MC)
-  MCEventHeader_2 mcevth_2[100];   //@< Event Header class (MC) for reflector > 0.6
+  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[100];         //@< number of photoelectrons in an event from showers
-  float inumphensb[100];             //@< number of photoelectrons in an event from nsb
+  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[100], phiCT[100];       //@< parameters of a given shower
+  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[100];              //@< impact parameter
+  float impactD[MAX_NUMBER_OF_CTS];  //@< impact parameter
   float l1, m1, n1;           //@< auxiliary variables
-  float l2, m2, n2;           //@< auxiliary variables
-  float num, den;      //@< auxiliary variables
-  float factorqe_NSB[100];          //@< factor on the NSB depending of QE
+  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[100];            //@< partial number of photons in a given run
-  int ntcph[100];                //@< total number of photons
+  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
@@ -375,10 +398,14 @@
   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[100][iMAXNUMPIX];  //@< diffuse NSB values for each pixel  
-                                     //@< derived from meanNSB
-  float nsbrate_phepns[iMAXNUMPIX][iNUMWAVEBANDS];    //@< non-diffuse nsb 
-                                                     //@< photoelectron rates 
-  float nsb_phepns[100][iMAXNUMPIX];                 //@< NSB values for each pixel
-  float nsb_phepns_rotated[100][iMAXNUMPIX];         //@< NSB values for each pixel after rotation
+  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
@@ -406,5 +433,5 @@
   int flagstoring = 0;
 
-  int ntrigger[100];           //@< number of triggers in the whole file
+  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
@@ -435,6 +462,4 @@
   Float_t telestheta = 0.0;
   Float_t telesphi = 0.0;
-  Float_t sofftheta = 0.0;
-  Float_t soffphi = 0.0;
   UInt_t corsika = 5200 ;
   Float_t maxpimpact = 0.0;
@@ -453,10 +478,30 @@
   //@'
 
-  qThreshold = new float *[100];
-  for(int i=0;i<100;i++)
+  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<100;ict++){      
+    for(int ict=0;ict<MAX_NUMBER_OF_CTS;ict++){      
       nsb_phepns[ict][i]=0;
       nsb_phepns_rotated[ict][i]=0;
@@ -465,5 +510,5 @@
       nsbrate_phepns[i][j]=0.0;    //@< Starlight rates initialised at 0.0
   }
-  for(int i=0;i<100;i++)
+  for(int i=0;i<MAX_NUMBER_OF_CTS;i++)
     ntcph[i]=0;
   /*!@'
@@ -527,13 +572,25 @@
 
   ct_Number=get_ct_number();
-  ct_Geometry=get_ct_geometry();
+  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;
-    if(ct_Geometry>=ict*10){
-      GeometryCamera[ict]=int(ct_Geometry/pow(10.0,ict))%10;
-    }
-    else
-      GeometryCamera[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);
   }
 
@@ -615,9 +672,11 @@
   get_secure_threhold(&riseDiskThres,&secureDiskThres);
 
-  get_FADC_properties( &FADC_response_ampl, &FADC_response_fwhm,
-		       &FADC_resp_ampl_out, &FADC_resp_fwhm_out);
+  get_FADC_properties
+    (&FADC_shape, &FADC_response_integ, &FADC_response_fwhm,
+     &FADC_shape_out, &FADC_resp_integ_out, &FADC_resp_fwhm_out);
+
   FADC_high2low=get_High_to_Low();
 
-  // FIXME --- tirgger properties may depend on the Camrea geometry
+  // FIXME --- trigger properties may depend on the Camera geometry
 
   get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm);
@@ -630,5 +689,5 @@
     (Trigger_loop_utop-Trigger_loop_ltop+1);
 
-  // Trigger loop operation is not implemented for Multi telscopes
+  // Trigger loop operation is not implemented for Multi telescopes
 
   if ( Trigger_Loop && ct_Number > 1 ){
@@ -651,5 +710,12 @@
 
   for(int ict=0;ict<ct_Number;ict++)
-    strcpy( inname_CT[ict], get_input_filename(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() );
@@ -662,6 +728,8 @@
   // get different parameters of the simulation
 
-  addElecNoise = add_elec_noise(&FADC_noise, &DIGITAL_noise, &Trigger_noise);
+  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);
 
 
@@ -681,4 +749,6 @@
 	"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
@@ -748,11 +818,13 @@
   
   log(SIGNATURE,
-      "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %3.2f(%s) %3.2f+%3.2f(%s) %s\n",
+      "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s %3.2f, %s %3.2f (%s), %3.2f (%s) + %3.2f (%s) %s\n",
       "Flags",
       "Data_From_STDIN",   ONoff(Data_From_STDIN),  
       "Write_All_Events",  ONoff(Write_All_Images), 
       "Rotate Starfield",  ONoff(Starfield_rotate),
-      "Electronic Noise",  Trigger_noise,"trigger", 
-      FADC_noise,DIGITAL_noise,"fadc",ONoff(addElecNoise)
+      "Electronic Noise",  "trigger (rel. noise): ", Trigger_noise, 
+      "FADC (ADC counts): ", FADC_noise_inner, "inner pixels",
+      FADC_noise_outer, "outer pixels", 
+      DIGITAL_noise, "added digital noise at FADC", ONoff(addElecNoise)
       );
 
@@ -779,5 +851,5 @@
 
   ntriggerloop= new int ** [(int)((Trigger_loop_uthres-Trigger_loop_lthres)
-			   /Trigger_loop_sthres)]; 
+			   /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];
@@ -821,5 +893,5 @@
 
   Int_t Lev0, Lev1; 
-  Int_t Lev0MT[100], Lev1MT[100]; 
+  Int_t Lev0MT[MAX_NUMBER_OF_CTS], Lev1MT[MAX_NUMBER_OF_CTS]; 
   
   fadcValues    =  new TArrayC(FADC_SLICES);
@@ -843,10 +915,10 @@
 
   for (int i=0; i<ct_Number;i++){
-    Trigger_CT[i] = new  MTrigger(TriggerPixels[i],
+    Trigger_CT[i] = new MTrigger(TriggerPixels[i],
 		    ((MGeomCam*)(camgeom.UncheckedAt(i))),
 		    Trigger_gate_length, 
 		    Trigger_overlaping_time,
 		    Trigger_response_ampl, 
-		    Trigger_response_fwhm);  //@< A instance of the Class MTrigger 
+		    Trigger_response_fwhm, i);  //@< A instance of the Class MTrigger 
     Trigger_CT[i]->SetSeed(UInt_t(i+get_seeds(0)));
   }
@@ -912,6 +984,8 @@
     Fadc_CT[i] = 
       new  MFadc(((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels(),
-		 FADC_response_ampl,FADC_response_fwhm,
-		 FADC_resp_ampl_out,FADC_resp_fwhm_out,
+		 FADC_shape,
+		 FADC_response_integ,FADC_response_fwhm,
+		 FADC_shape_out,
+		 FADC_resp_integ_out,FADC_resp_fwhm_out,
 		 get_trig_delay()
 	  ) ; //@< A instance of the Class MFadc
@@ -940,8 +1014,19 @@
     
     // Generate database for the Fadc electronic noise
-    if (addElecNoise){        
-      Fadc_CT[ict]->SetElecNoise(FADC_noise);
-      Fadc_CT[ict]->SetDigitalNoise(DIGITAL_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);
   }
 
@@ -955,5 +1040,10 @@
   MRawRunHeader *RunHeader= new MRawRunHeader();
   MMcRunHeader  *McRunHeader = new MMcRunHeader();
-  MMcCorsikaRunHeader  *McCorsikaRunHeader = new MMcCorsikaRunHeader();
+  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;
@@ -973,4 +1063,5 @@
       EvtHeader[i] = new MRawEvtHeader();
     }
+
   }
 
@@ -1015,26 +1106,27 @@
   char help[4];  
 
-  HeaderTree.Branch("MRawRunHeader","MRawRunHeader",
+  HeaderTree.Branch("MRawRunHeader.","MRawRunHeader",
   		    &RunHeader);
 
-  HeaderTree.Branch("MMcRunHeader","MMcRunHeader",
+  HeaderTree.Branch("MMcRunHeader.","MMcRunHeader",
   		    &McRunHeader);
 
-  HeaderTree.Branch("MMcCorsikaRunHeader","MMcCorsikaRunHeader",
-  		    &McCorsikaRunHeader);
+  HeaderTree.Branch("MMcCorsikaRunHeader.","MMcCorsikaRunHeader",
+		    &McCorsikaRunHeader);
 
   if(!Trigger_Loop && Write_McTrig && ct_Number==1){
     
-    HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader", 
+    HeaderTree.Branch("MMcTrigHeader.","MMcTrigHeader", 
 		      &HeaderTrig[0]);
   }
+
   if (ct_Number==1){
-    HeaderTree.Branch("MGeomCam",GeometryName[0],
-		      &camdummy[0]);
-    HeaderTree.Branch("MMcConfigRunHeader","MMcConfigRunHeader",
+    //    HeaderTree.Branch("MGeomCam.", "MGeomCam", &camdummy[0]);
+    HeaderTree.Branch("MGeomCam.", GeometryName[0], &camdummy[0]);
+    HeaderTree.Branch("MMcConfigRunHeader.","MMcConfigRunHeader",
 		      &McConfigRunHeader[0]);
   }
   else{
-    char branchname[20];
+    char branchname[256];
     for (int ict=0; ict<ct_Number;ict++){
       sprintf(help,"%i",ict+1);
@@ -1042,5 +1134,5 @@
       strcat (branchname, & help[0]);
       strcat (branchname, ".");
-      HeaderTree.Branch(branchname,GeometryName[ict], 
+      HeaderTree.Branch(branchname, GeometryName[ict], 
 			&camdummy[ict]);
     }
@@ -1058,5 +1150,5 @@
   if ((Trigger_Loop || ct_Number>1) && Write_McTrig){
     ibr=0;
-    for(char branchname[20];ibr<numberBranches;ibr++){
+    for(char branchname[256];ibr<numberBranches;ibr++){
       
       sprintf(help,"%i",ibr+1);
@@ -1071,10 +1163,10 @@
   if(!Trigger_Loop && Write_McFADC && ct_Number==1){
     
-    HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader", 
+    HeaderTree.Branch("MMcFadcHeader.","MMcFadcHeader", 
 		      &HeaderFadc[0]);
   }
   if ((Trigger_Loop || ct_Number>1) && Write_McFADC){
     ibr=0;
-    for(char branchname[20];ibr<numberBranches;ibr++){
+    for(char branchname[256];ibr<numberBranches;ibr++){
       
       sprintf(help,"%i",ibr+1);
@@ -1082,4 +1174,5 @@
       strcat (branchname, & help[0]);
       strcat (branchname, ".");
+
       HeaderTree.Branch(branchname,"MMcFadcHeader", 
 			&HeaderFadc[ibr]);
@@ -1094,5 +1187,5 @@
   RunHeader->SetRunType(256);
   RunHeader->SetRunNumber(0);
-  RunHeader->SetNumSamples(0,FADC_SLICES);
+  RunHeader->SetNumSamples(FADC_SLICES, FADC_SLICES);
   RunHeader->SetNumCrates(1);
   RunHeader->SetNumPixInCrate(ct_NPixels);
@@ -1149,21 +1242,21 @@
   
   //  Fill branches for MMcFadcHeader
-
-  for(int i=0;i<ct_NPixels;i++){
-      fadc_elecnoise[i]=FADC_noise;
-      fadc_diginoise[i]=DIGITAL_noise;
-  }
-
   Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
 
   if(!Trigger_Loop && Write_McFADC && ct_Number==1){
 
-    HeaderFadc[0]->SetShape(0.0);
-    HeaderFadc[0]->SetAmplitud(FADC_response_ampl*2.35/
-			       sqrt(2*M_PI*FADC_response_fwhm
-				    *FADC_response_fwhm),
-			       FADC_resp_ampl_out*2.35/
-			       sqrt(2*M_PI*FADC_resp_fwhm_out
-				    *FADC_resp_fwhm_out));
+      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(FADC_shape);
+    HeaderFadc[0]->SetShapeOuter(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);
@@ -1173,14 +1266,22 @@
 				((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
   }
+
   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(0.0);
-      HeaderFadc[i]->SetAmplitud(FADC_response_ampl*2.35/
-				 sqrt(2*M_PI*FADC_response_fwhm
-				      *FADC_response_fwhm),
-				 FADC_resp_ampl_out*2.35/
-				 sqrt(2*M_PI*FADC_resp_fwhm_out
-				      *FADC_resp_fwhm_out));
+      HeaderFadc[i]->SetShape(FADC_shape);
+      HeaderFadc[i]->SetShapeOuter(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);
@@ -1189,6 +1290,17 @@
 				  ((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
     }
-  }  
+  }
+  
   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){
@@ -1196,11 +1308,8 @@
 	for(itopocount=0;itopocount<=Trigger_loop_utop-Trigger_loop_ltop;itopocount++){
 	  Fadc_CT[0]->GetPedestals(&fadc_pedestals[0]);
-	  HeaderFadc[iconcount]->SetShape(0.0);
-	  HeaderFadc[iconcount]->SetAmplitud(FADC_response_ampl*2.35/
-					     sqrt(2*M_PI*FADC_response_fwhm
-						  *FADC_response_fwhm),
-					     FADC_resp_ampl_out*2.35/
-					     sqrt(2*M_PI*FADC_resp_fwhm_out
-						  *FADC_resp_fwhm_out));
+	  HeaderFadc[iconcount]->SetShape(FADC_shape);
+	  HeaderFadc[iconcount]->SetShapeOuter(FADC_shape_out);
+	  HeaderFadc[iconcount]->SetAmplitud(FADC_response_integ,
+					     FADC_resp_integ_out);
 	  HeaderFadc[iconcount]->SetFwhm(FADC_response_fwhm,
 					 FADC_resp_fwhm_out);
@@ -1222,5 +1331,4 @@
   } 
 
-
   //      create a Tree for the Event data stream 
   TTree EvtTree("Events","Normal Triggered Events");
@@ -1228,5 +1336,5 @@
   if (Write_McEvt && ct_Number==1){
 
-    EvtTree.Branch("MMcEvt","MMcEvt", 
+    EvtTree.Branch("MMcEvt.","MMcEvt", 
 		   &McEvt[0]);
   }
@@ -1234,5 +1342,5 @@
   if (Write_McEvt && ct_Number!=1){
     ibr=0;
-    for(char branchname[10];ibr<numberBranches;ibr++){
+    for(char branchname[256];ibr<numberBranches;ibr++){
       
       sprintf(help,"%i",ibr+1);
@@ -1248,11 +1356,11 @@
     
     if (Write_RawEvt){
-      EvtTree.Branch("MRawEvtHeader","MRawEvtHeader", 
+      EvtTree.Branch("MRawEvtHeader.","MRawEvtHeader", 
 		     &EvtHeader[0]);
-      EvtTree.Branch("MRawEvtData","MRawEvtData", 
+      EvtTree.Branch("MRawEvtData.","MRawEvtData", 
 		     &EvtData[0]);
     }
     if (Write_McTrig){
-      EvtTree.Branch("MMcTrig","MMcTrig", 
+      EvtTree.Branch("MMcTrig.","MMcTrig", 
 		     &McTrig[0]);
     }    
@@ -1262,5 +1370,5 @@
     if (Write_McTrig){
       ibr=0;
-      for(char branchname[10];ibr<numberBranches;ibr++){
+      for(char branchname[256];ibr<numberBranches;ibr++){
       
 	sprintf(help,"%i",ibr+1);
@@ -1276,5 +1384,5 @@
   if ((Trigger_Loop || ct_Number>1) && Write_RawEvt){
     ibr=0;
-    for(char branchname[15];ibr<numberBranches;ibr++){
+    for(char branchname[256];ibr<numberBranches;ibr++){
       
       sprintf(help,"%i",ibr+1);
@@ -1286,5 +1394,5 @@
     }
     ibr=0;
-    for(char branchname[15];ibr<numberBranches;ibr++){
+    for(char branchname[256];ibr<numberBranches;ibr++){
       
       sprintf(help,"%i",ibr+1);
@@ -1315,10 +1423,9 @@
   }
   
-
   // prepare the NSB simulation
 
   //  Instance of the Mlons class
-  MLons lons(Trigger_response_ampl, Trigger_response_fwhm,
-	     FADC_response_ampl,FADC_response_fwhm);
+  MLons lons(0.0, Trigger_response_ampl, Trigger_response_fwhm,
+	     float(FADC_shape),FADC_response_integ,FADC_response_fwhm);
 
   lons.SetSeed(Int_t(get_seeds(1)/get_seeds(0))+1);
@@ -1327,6 +1434,6 @@
 
   //  Instance of the Mlons class
-  MLons lons_outer(Trigger_response_ampl, Trigger_response_fwhm,
-		   FADC_resp_ampl_out,FADC_resp_fwhm_out);
+  MLons lons_outer(0.0, Trigger_response_ampl, Trigger_response_fwhm,
+		   float(FADC_shape_out),FADC_resp_integ_out,FADC_resp_fwhm_out);
 
   lons_outer.SetSeed(Int_t(get_seeds(1)/get_seeds(0))+2);
@@ -1340,5 +1447,8 @@
     //
     
-    // FIXME --- star NSB different for each camera?
+    //
+    // FIXME! --- star NSB different for each camera?
+    // Then we will have to use mirror_frac[ict]
+    //
     log(SIGNATURE,"Produce NSB rates from Star Field");
 
@@ -1346,5 +1456,16 @@
 			  ((MGeomCam*)(camgeom.UncheckedAt(0))),
 			  nsbrate_phepns,
-			  0);
+			  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){
@@ -1358,5 +1479,6 @@
     for(int ict=0;ict<ct_Number;ict++){      
       
-      // First we set the factor due to the mirror size respect to the MAGIC mirror.
+      // First we set the factor due to the mirror size with respect to the normal 
+      // MAGIC mirror geometry.
       
       switch(GeometryCamera[ict]){
@@ -1379,14 +1501,31 @@
       }
 
-      factorNSB_ct=factorNSB_ct/
+      factorNSB_ct = factorNSB_ct /
 	((*((MGeomCam*)(camgeom.UncheckedAt(ict)))).GetCameraDist()*1000*
 	 (*((MGeomCam*)(camgeom.UncheckedAt(ict)))).GetCameraDist()*1000*
-	 PIXEL_SIZE*PIXEL_SIZE);
+	 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();
-	diffnsb_phepns[ict][ui] = (Int_t(meanNSB*factorNSB_ct*100*size*size+0.5))/(100.0)*factorqe_NSB[ict];
+	  (*((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];
       }
     }
@@ -1396,5 +1535,5 @@
     for(int ict=0;ict<ct_Number;ict++){      
       for(j=0;j<iNUMWAVEBANDS;j++){
-	// calculate the effect of the atmospheric extinction 
+	// calculate the effect of the atmospheric extinction (for stars!)
 	
 	zenfactor = pow(10., -0.4 * ext[j] );
@@ -1411,12 +1550,38 @@
   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   //
-  //  Now 500 empty event with the conditin in which the camera is run
-  //  are simulated. IN this way one gets an estimation of the 
-  //  sigma for the pedestal of each FADC channel.
-  //  This sigma computtaion is done assuming any noise taht effects 
+  //  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)
+  //  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++){
@@ -1427,8 +1592,8 @@
       fadc_sigma_low[ui]=0.0;
     }
-    for(int ie=0;ie<500;ie++){
+    for(int ie=0; ie < empty_events; ie++){
       Fadc_CT[ict]->Reset() ; 
       if (addElecNoise){	
-	Fadc_CT[ict]->ElecNoise(FADC_noise) ;
+	Fadc_CT[ict]->ElecNoise() ;
       }
       if(simulateNSB){
@@ -1455,11 +1620,36 @@
       }
       Fadc_CT[ict]->Pedestals();
+
       for(UInt_t ui=0;
 	  ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
-	  ui++){
-	fadc_sigma[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,1))/500.0;
-	fadc_sigma_low[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,0))/500.0;
+	  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)))
@@ -1560,5 +1750,9 @@
 	}
 
-	fread((char*)&mcrunh,(SIZE_OF_MCRUNHEADER)*sizeof(float), 1, inputfile[ict] );
+	//	fread((char*)&mcrunh,(SIZE_OF_MCRUNHEADER)*sizeof(float), 1, inputfile[ict] );
+	//  AM, changed reading method, 02/2004:
+
+	mcrunh.read(inputfile[ict]);
+
       }
       
@@ -1598,41 +1792,34 @@
 	if (reflector_file_version<6){
 	  for(int ict=0;ict<ct_Number;ict++)
-	    fread( (char*)&mcevth[ict], mcevth[ict].mysize(),
-		   1, inputfile[ict] );
+	    {
+	      mcevth[ict].read(inputfile[ict]);
+	    }
+
 	}
 	else{
 	  for(int ict=0;ict<ct_Number;ict++)
-	    fread( (char*)&mcevth_2[ict], mcevth_2[ict].mysize(),
-		   1, inputfile[ict] );
+	    {
+	      mcevth_2[ict].read(inputfile[ict]);
+	    }
 	}
 
-	// calculate impact parameter (shortest distance betwee the original
-	// trajectory of the primary (assumed shower-axis) and the
-	// direction where the telescope points to
+	//
+	// 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.
 	// 
-	// we use the following equation, given that the shower core position
-	// is (x1,y1,z1)=(x,y,0),the  trajectory is given by (l1,m1,n1),
-	// and the telescope position and orientation are (x2,y2,z2)=(0,0,0)
-	// and (l2,m2,n2)
-	//
-	//               |                     |
-	//               | x1-x2  y1-y2  z1-z2 |
-	//               |                     |
-	//             + |   l1     m1     n1  |
-	//             - |                     |
-	//               |   l2     m2     n2  |
-	//               |                     |
-	// dist = ------------------------------------        ( > 0 )
-	//        [ |l1 m1|2   |m1 n1|2   |n1 l1|2 ]1/2
-	//        [ |     |  + |     |  + |     |  ]
-	//        [ |l2 m2|    |m2 n2|    |n2 l2|  ]
-	//
-	// playing a little bit, we get this reduced for in our case:
-	//
-	//
-	// dist = (- m2 n1 x + m1 n2 x + l2 n1 y - l1 n2 y - l2 m1 z + l1 m2 z) /
-	//         [(l2^2 (m1^2 + n1^2) + (m2 n1 - m1 n2)^2 - 
-	//          2 l1 l2 (m1 m2 + n1 n2) + l1^2 (m2^2 + n2^2) ] ^(1/2)
-	
+
+
 	// read the direction of the incoming shower
 	// It is done only for one telescope since it is suposed
@@ -1646,4 +1833,5 @@
 	  phishw = mcevth_2[0].get_phi();
 	}
+	Zenith=thetashw; Azimutal=phishw;
 	
 	// calculate vector for shower
@@ -1653,79 +1841,31 @@
 	n1 = cos(thetashw);
 
-	if (reflector_file_version<6){
-
-	  mcevth[0].get_core(&coreX, &coreY);	
-	  
-	  for(int ict=0;ict<ct_Number;ict++){
-	    // read the deviation of the telescope with respect to the shower
-	    mcevth[ict].get_deviations ( &thetaCT[ict], &phiCT[ict] );
-
-	    if ( (thetaCT[ict] == 0.) && (phiCT[ict] == 0.) ) {
-	  
-	      // CT was looking to the source (both lines are parallel)
-	      // therefore, we calculate the impact parameter as the distance 
-	      // between the CT axis and the core position
-
-	      impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
-	      thetaCT[ict] += thetashw;
-	      phiCT[ict]   += phishw;
-	    } else {
-	    
-	      // the shower comes off-axis
-	    
-	      // obtain with this the final direction of the CT
-	      
-	      thetaCT[ict] += thetashw;
-	      phiCT[ict]   += phishw;
-	      
-	      // calculate vector for telescope
-	    
-	      l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
-	      m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
-	      n2 = cos(thetaCT[ict]);
-	      
-	      num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
-	      den = (SQR(l1*m2 - l2*m1) + 
-		     SQR(m1*n2 - m2*n1) + 
-		     SQR(n1*l2 - n2*l1));
-	      den = sqrt(den);
-	      impactD[ict] = fabs(num)/den;
-	    }
-	  }
-	}
-	else{
-	  mcevth_2[0].get_core(&coreX, &coreY);	
-	  for(int ict=0;ict<ct_Number;ict++){
-	    thetaCT[ict]=mcevth_2[ict].get_theta_CT();
-	    phiCT[ict]=mcevth_2[ict].get_phi_CT();
-	    if ( (thetaCT[ict] == thetashw) && (phiCT[ict] == phishw) ) {
-	  
-	      // CT was looking to the source (both lines are parallel)
-	      // therefore, we calculate the impact parameter as the distance 
-	      // between the CT axis and the core position
-
-	      impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
-	    
-	    } else {
-	    
-	      // the shower comes off-axis
-	    
-	      // calculate vector for telescope
-	    
-	      l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
-	      m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
-	      n2 = cos(thetaCT[ict]);
-	    
-	      num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
-	      den = (SQR(l1*m2 - l2*m1) + 
-		     SQR(m1*n2 - m2*n1) + 
-		     SQR(n1*l2 - n2*l1));
-	      den = sqrt(den);
-	    
-	      impactD[ict] = fabs(num)/den;
-	    
-	    }
-	  }
-	}
+	//
+	// 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
@@ -1772,5 +1912,6 @@
 			    &arrtmin_ns, // will be changed by the function!
 			    &arrtmax_ns, // will be changed by the function!
-			    ict);
+			    ict,
+			    mirror_frac[ict]);
 
 	  inumphe=(inumphe<inumphe_CT[ict])?inumphe_CT[ict]:inumphe;
@@ -1908,5 +2049,5 @@
 	    Trigger_CT[ict]->ElecNoise(Trigger_noise) ;
 	    
-	    Fadc_CT[ict]->ElecNoise(FADC_noise) ;
+	    Fadc_CT[ict]->ElecNoise() ;
 	  }
 	}
@@ -2079,6 +2220,6 @@
 			     mcevth[0].get_phi(), 
 			     mcevth[0].get_core(),
-			     mcevth[0].get_coreX(),
-			     mcevth[0].get_coreY(),
+			     coreX,
+			     coreY,
 			     impactD[0],
 			     phiCT[0],
@@ -2116,6 +2257,6 @@
 			     mcevth_2[0].get_phi(), 
 			     mcevth_2[0].get_core(),
-			     mcevth_2[0].get_coreX(),
-			     mcevth_2[0].get_coreY(),
+			     coreX,
+			     coreY,
 			     impactD[0],
 			     mcevth_2[0].get_phi_CT(),
@@ -2314,4 +2455,5 @@
 	      //
 	    
+
 	      for (int ict=0;ict<ct_Number;ict++){
 		Float_t ftime, ltime;
@@ -2327,6 +2469,6 @@
 				  mcevth[ict].get_phi(), 
 				  mcevth[ict].get_core(),
-				  mcevth[ict].get_coreX(),
-				  mcevth[ict].get_coreY(),
+				  coreX,
+				  coreY,
 				  impactD[ict],
 				  phiCT[ict],
@@ -2364,6 +2506,6 @@
 				  mcevth_2[ict].get_phi(), 
 				  mcevth_2[ict].get_core(),
-				  mcevth_2[ict].get_coreX(),
-				  mcevth_2[ict].get_coreY(),
+				  coreX,
+				  coreY,
 				  impactD[ict],
 				  mcevth_2[ict].get_phi_CT(),
@@ -2470,11 +2612,14 @@
 	  for(int ict=0;ict<ct_Number;ict++){
 	    read_ascii(inputfile[ict], McConfigRunHeader[ict]);
-
-	    if (get_sigma_xy_cm_spot() > 0.)
+	    McConfigRunHeader[ict]->SetMissPointingX(missP_x);
+	    McConfigRunHeader[ict]->SetMissPointingY(missP_y);
+	    if ( Spotsigma > 0.)
 	    {
-		Float_t spotsigma = McConfigRunHeader[ict]->GetPointSpread();
-		Float_t newsigma = sqrt(spotsigma * spotsigma + 
-		    get_sigma_xy_cm_spot() * get_sigma_xy_cm_spot());
-		McConfigRunHeader[ict]->SetPointSpread(newsigma);
+		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);
 	    }
 	  }
@@ -2534,5 +2679,4 @@
   }
 
-  get_source_off(&sofftheta,&soffphi);
   if(!Trigger_Loop) icontrigger=0;
   time (&ltime);
@@ -2564,6 +2708,4 @@
 		    telestheta,
 		    telesphi,
-		    sofftheta,
-		    soffphi,
 		    shthetamax,
 		    shthetamin,
@@ -2605,6 +2747,4 @@
 		    telestheta,
 		    telesphi,
-		    sofftheta,
-		    soffphi,
 		    shthetamax,
 		    shthetamin,
@@ -2623,5 +2763,5 @@
   // Fill some missing values for MRawRunHeader
 
-  RunHeader->SetRunNumber(rnum);
+  RunHeader->SetRunNumber((UInt_t)rnum);
   RunHeader->SetNumEvents(nstoredevents);
 
@@ -2688,10 +2828,15 @@
     }
     
-    // Store Light Guides factors in the output file
-    TArrayF theta_lg;
-    TArrayF factor_lg;
-    theta_lg.Set(pointsWC,WC[0]);
-    factor_lg.Set(pointsWC,WC[1]);
-    McConfigRunHeader[ict]->SetLightGuides(theta_lg,factor_lg);
+    // 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);
   
   }
@@ -2980,8 +3125,8 @@
       break;
 
-    if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
+    if ( (i < ct_NPixels) && (i > -1) &&
          ((j-1) < pointsQE[ict])   && ((j-1) > -1) ) {
-      QE[ict][i-1][0][j-1] = QElambda[j-1];
-      QE[ict][i-1][1][j-1] = qe;
+      QE[ict][i][0][j-1] = QElambda[j-1];
+      QE[ict][i][1][j-1] = qe;
     }
 
@@ -3036,5 +3181,5 @@
 
   //------------------------------------------------------------
-  // Read Light Guides data
+  // Read Light Collection data
 
   // try to open the file
@@ -3055,6 +3200,8 @@
   // get line from the file
 
-  wcfile.getline(line, LINE_MAX_LENGTH);
-  
+  do
+    wcfile.getline(line, LINE_MAX_LENGTH);
+  while (line[0] == '#');
+
   // get the number of datapoints 
     
@@ -3071,5 +3218,29 @@
     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
 
@@ -3548,5 +3719,7 @@
 // dist_r_P
 //
-// distance straight line r - point 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
 //------------------------------------------------------------
 
@@ -3559,7 +3732,5 @@
           sqrt((SQR((a-x)*m-(b-y)*l) +
                 SQR((b-y)*n-(c-z)*m) +
-                SQR((c-z)*l-(a-x)*n))/
-               (SQR(l)+SQR(m)+SQR(n))
-               )
+                SQR((c-z)*l-(a-x)*n)))
           );
 }
@@ -3636,5 +3807,5 @@
 int size_rotated(
     float *rotated,
-    float nsb[iMAXNUMPIX],
+    float *nsb,
     float rho)
 {
@@ -3737,15 +3908,16 @@
 		  class MFadc *fadc,
 		  int *itotnphe, // total number of produced photoelectrons
-		  float nphe[iMAXNUMPIX], // number of photoelectrons in each pixel
+		  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.
+		  int ict,      // Telescope that is being analised to get the right QE.
+		  float mirror_fraction // Fraction of total mirror really present
                    ){
 
   static uint i;
   static int k, ipixnum;
-  static float cx, cy, wl, qe, t;
-  static float cu, cv, cw;
+  static float cx, cy, wl, qe;
+  static float cw;
   static MCCphoton photon;
   static float **qept;
@@ -3753,4 +3925,6 @@
   static float radius_mm;
 
+  Float_t t;
+
   // reset variables
 
@@ -3761,4 +3935,5 @@
   }
 
+  *itotnphe = 0;
   *incph = 0;
 
@@ -3768,5 +3943,6 @@
   random.SetSeed((Int_t)(get_seeds(0)*get_seeds(1)));
 
-  float spotsigma = get_sigma_xy_cm_spot();
+  float C1, C2, C3, rho;
+
 
   //- - - - - - - - - - - - - - - - - - - - - - - - - 
@@ -3779,28 +3955,35 @@
   
   // read the photons data
-  
-  
+
   // loop over the photons
   
   while (1) {
 
-    fread ( flag, SIZE_OF_FLAGS, 1, sp );
+    photon.read(sp);
+    photon.get_flag(flag);
 
     if (isA( flag, FLAG_END_OF_EVENT ))
-        break;
-	  
-    memcpy( (char*)&photon, flag, SIZE_OF_FLAGS );
-    
-    fread( ((char*)&photon)+SIZE_OF_FLAGS, photon.mysize()-SIZE_OF_FLAGS, 1, sp );
-	  
+      {
+	fseek (sp, SIZE_OF_FLAGS-photon.mysize(), SEEK_CUR);
+	break;  
+      }
 
     //  Check if photon is inside trigger time range
 
-    t = photon.get_t() ; 
+    t = photon.get_t() ;
 
     if (t-*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
@@ -3810,10 +3993,24 @@
     cy = photon.get_y(); 
 
-
-    if (spotsigma > 0.)
+    if (Spotsigma > 0.)
     {
-	cx += random.Gaus(0.,spotsigma);
-	cy += random.Gaus(0.,spotsigma);
-    }
+	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
@@ -3849,5 +4046,5 @@
     
     // set pointer to the QE table of the relevant pixel
-    	  
+
     qept = (float **)QE[ict][ipixnum];
 	  
@@ -3870,20 +4067,32 @@
     qe = lin_interpol(qept[0][k-1], qept[1][k-1], qept[0][k], qept[1][k], wl) / 100.0;
 
-    // Apply incient angular correction due to Light Guides
-
-    cu=photon.get_u();
-    cv=photon.get_v();
-    cw=90.0-acos(sqrt(1-cu*cu-cv*cv))*3.14159/180.0;
-
-    // find data point in the QE table (-> k)
-
-    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);
+    //
+    // Apply incient 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
@@ -3902,5 +4111,6 @@
 	  
     nphe[ipixnum] += 1.0;
-	  
+
+
     // store the new photoelectron
 
@@ -3929,9 +4139,10 @@
 int produce_nsbrates( char *iname, // the starfield input file name
 		      MGeomCam *camgeom, // camera layout
-		      float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS], // the product of this function:
-		                               // the NSB rates in phe/ns for each pixel 
-		      int ict
-		      ){
-
+		      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
@@ -3959,4 +4170,5 @@
   char flag_new[4]; 
 
+
   // open input file 
   
@@ -3965,5 +4177,5 @@
   infile = fopen( iname, "r" );
 
-  // check if the satrfield input file exists
+  // check if the starfield input file exists
 
   if ( infile == NULL ) {
@@ -3975,5 +4187,5 @@
     return (0);
   }
-  
+
   // get signature, and check it
   
@@ -3982,11 +4194,24 @@
   }
 
-  // Instance of MTrigger and MFadc
-
-  MTrigger trigger(camgeom->GetNumPixels(),camgeom,Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
-  MFadc flashadc(camgeom->GetNumPixels());
-
-  // initialize flag
-  
+  // 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());
+
+  // initialize flag  
   strcpy( cflag, "                                        \0" );
 
@@ -4038,8 +4263,10 @@
 	
 	if (reflector_file_version<6)
-	  fread( (char*)&evth, evth.mysize(), 1, infile );
+	  //      fread( (char*)&evth, evth.mysize(), 1, infile );
+	  evth.read(infile);
 	else
-	  fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
-	  
+	  //	  fread( (char*)&evth_2, evth_2.mysize(), 1, infile );
+	  evth_2.read(infile);
+
 	if (reflector_file_version<6)
 	  integtime_ns = evth.get_energy();
@@ -4081,6 +4308,6 @@
 			      wl_nm[i],
 			      wl_nm[i+1],
-			      &trigger,  // this is a dummy here
-			      &flashadc, // this is a dummy here
+			      trigger,  // this is a dummy here
+			      flashadc, // this is a dummy here
 			      &itnphe,
 			      nphe, // we want this!
@@ -4088,6 +4315,6 @@
 			      &tmin_ns,
 			      &tmax_ns,
-			      0
-			      ); // photons from starfield
+			      0,
+			      mirror_fraction); // photons from starfield
 
 	    if( k != 0 ){ // non-zero returnvalue means error
@@ -4124,125 +4351,4 @@
 
   return(0);
-
-}
-
-
-//------------------------------------------------------------
-// @name produce_nsb_phes
-//                                     
-// @desc produce the photoelectrons from the nsbrates
-//
-// @date Thu Feb 17 17:10:40 CET 2000
-// @function @code 
-//------------------------------------------------------------
-
-int produce_nsb_phes( float *atmin_ns,
-		      float *atmax_ns,
-		      float theta_rad,
-		      struct camera *cam,
-		      float nsbr_phepns[iMAXNUMPIX][iNUMWAVEBANDS],
-		      float difnsbr_phepns[iMAXNUMPIX],
-		      float extinction[iNUMWAVEBANDS],
-		      float fnpx[iMAXNUMPIX],
-		      MTrigger *trigger,
-		      MFadc *fadc,
-		      int *inphe,
-		      float base_mv[iMAXNUMPIX]){
-
-  float simtime_ns; // NSB simulation time
-  int i, j, k, ii; 
-  float zenfactor;  //  correction factor calculated from the extinction
-  int inumnsbphe;   //  number of photoelectrons caused by NSB
-
-  float t;
-  TRandom random;             // Random numbers generator
-
-  random.SetSeed(get_seeds(1));
-
-  ii = *inphe; // avoid dereferencing
-		
-  // check if the arrival times are set; if not generate them
-
-  if(*atmin_ns <SIMTIMEOFFSET_NS || *atmin_ns > *atmax_ns){
-    *atmin_ns = 0.;
-    *atmax_ns = simtime_ns = TOTAL_TRIGGER_TIME;
-
-  }
-  else{ // extend the event time window by the given offsets
-
-    *atmin_ns = *atmin_ns - SIMTIMEOFFSET_NS;
-    *atmax_ns = *atmax_ns + SIMTIMEOFFSET_NS;
-
-    simtime_ns = *atmax_ns - *atmin_ns;
-
-    // make sure the simulated time is long enough for the FADC 
-    // simulation and not too long
-
-    if(simtime_ns< TOTAL_TRIGGER_TIME){
-      *atmin_ns = *atmin_ns -(TOTAL_TRIGGER_TIME-simtime_ns)/2;
-      *atmax_ns = *atmin_ns + TOTAL_TRIGGER_TIME;
-      simtime_ns = TOTAL_TRIGGER_TIME;
-    }
-
-    if(simtime_ns> TOTAL_TRIGGER_TIME){
-      *atmax_ns =*atmin_ns + TOTAL_TRIGGER_TIME;
-      simtime_ns = TOTAL_TRIGGER_TIME;
-    }
- 
-  }
-
-  // initialize baselines
-
-  for(i=0; i<cam->inumpixels; i++){
-    base_mv[i] = 0.;
-  }
-
-  // calculate baselines and generate phes
-
-  for(i=0; i<iNUMWAVEBANDS; i++){ // loop over the wavebands
-    
-    // calculate the effect of the atmospheric extinction 
-    
-    zenfactor = pow(10., -0.4 * extinction[i]/cos(theta_rad) );
-    
-    for(j=0; j<cam->inumpixels; j++){ // loop over the pixels
-      
-      inumnsbphe = (int) ((zenfactor * nsbr_phepns[j][i] + difnsbr_phepns[j]/iNUMWAVEBANDS)  
-			  * simtime_ns );
-
-      base_mv[j] += inumnsbphe;
-
-      // randomize
-
-      if (inumnsbphe>0.0){
-	inumnsbphe = random.Poisson (inumnsbphe );
-      }
-
-      // create the photoelectrons
-      
-      for(k=0; k < inumnsbphe; k++){
-	
-	t=(RandomNumber * simtime_ns);
-
-	(*fadc).Fill(j,t ,(*trigger).FillNSB(j,t));
-	
-	ii++; // increment total number of photoelectons
-	
-	fnpx[j] += 1.; // increment number of photoelectrons in this pixel
-
-      }
-      
-    } // end for(j=0 ...
-  } // end for(i=0 ...
-
-  // finish baseline calculation
-
-  for(i=0; i<cam->inumpixels; i++){
-    base_mv[i] *= RESPONSE_FWHM * RESPONSE_AMPLITUDE / simtime_ns;
-  }
-
-  *inphe = ii; // update the pointer
-
-  return(0);
 }
 
@@ -4257,4 +4363,28 @@
 //
 // $Log: not supported by cvs2svn $
+//
+// Revision 1.70  2004/09/15 A. 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]
@@ -4265,5 +4395,5 @@
 // has been introduced.
 //
-// The possibilty to unlarge the point spread function has been introduced
+// The possibilty to enlarge the point spread function has been introduced
 // in order to be able to simualte the current data.
 //
