Changeset 5072 for trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx

Timestamp:

09/16/04 16:16:34 (20 years ago)

Author:

moralejo

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx (modified) (99 diffs)

trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx

@@ -21 +21 @@
 //
 // $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 +84 @@
 #include "MGeomCamCT1Daniel.h"
 #include "MGeomPix.h"
+#include "MGeomCorsikaCT.h"
 
 /*!@" 
@@ -151 +152 @@
 static int ct_Number;        
 
-//@: Camera geometries
-static int ct_Geometry;        
-
 //@: list of showers to be skipped
 static int *Skip;
@@ -200 +198 @@
 
 //@: 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 +211 @@
 //@: 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 +223 @@
 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 +275 @@
 
 //@: 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 +330 @@
   //@'
 
-  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 +345 @@
 
   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 +398 @@
   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 +433 @@
   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 +462 @@
   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 +478 @@
   //@'
 
-  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 +510 @@
       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 +572 @@
 
   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 +672 @@
   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 +689 @@
     (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 +710 @@
 
   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 +728 @@
   // 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 +749 @@
         "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 +818 @@
   
   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 +851 @@
 
   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 +893 @@
 
   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 +915 @@
 
   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 +984 @@
     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 +1014 @@
     
     // 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 +1040 @@
   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 +1063 @@
       EvtHeader[i] = new MRawEvtHeader();
     }
+
   }
 
@@ -1015 +1106 @@
   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 +1134 @@
       strcat (branchname, & help[0]);
       strcat (branchname, ".");
-      HeaderTree.Branch(branchname,GeometryName[ict], 
+      HeaderTree.Branch(branchname, GeometryName[ict], 
                         &camdummy[ict]);
     }
@@ -1058 +1150 @@
   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 +1163 @@
   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 +1174 @@
       strcat (branchname, & help[0]);
       strcat (branchname, ".");
+
       HeaderTree.Branch(branchname,"MMcFadcHeader", 
                         &HeaderFadc[ibr]);
@@ -1094 +1187 @@
   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 +1242 @@
   
   //  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 +1266 @@
                                 ((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 +1290 @@
                                   ((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 +1308 @@
         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 +1331 @@
   } 
 
-
   //      create a Tree for the Event data stream 
   TTree EvtTree("Events","Normal Triggered Events");
@@ -1228 +1336 @@
   if (Write_McEvt && ct_Number==1){
 
-    EvtTree.Branch("MMcEvt","MMcEvt", 
+    EvtTree.Branch("MMcEvt.","MMcEvt", 
                    &McEvt[0]);
   }
@@ -1234 +1342 @@
   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 +1356 @@
     
     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 +1370 @@
     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 +1384 @@
   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 +1394 @@
     }
     ibr=0;
-    for(char branchname[15];ibr<numberBranches;ibr++){
+    for(char branchname[256];ibr<numberBranches;ibr++){
       
       sprintf(help,"%i",ibr+1);
@@ -1315 +1423 @@
   }
   
-
   // 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 +1434 @@
 
   //  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 +1447 @@
     //
     
-    // 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 +1456 @@
                           ((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 +1479 @@
     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 +1501 @@
       }
 
-      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 +1535 @@
     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 +1550 @@
   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   //
-  //  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 +1592 @@
       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 +1620 @@
       }
       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 +1750 @@
         }
 
-        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 +1792 @@
         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 +1833 @@
           phishw = mcevth_2[0].get_phi();
         }
+        Zenith=thetashw; Azimutal=phishw;
         
         // calculate vector for shower
@@ -1653 +1841 @@
         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 +1912 @@
                             &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 +2049 @@
             Trigger_CT[ict]->ElecNoise(Trigger_noise) ;
             
-            Fadc_CT[ict]->ElecNoise(FADC_noise) ;
+            Fadc_CT[ict]->ElecNoise() ;
           }
         }
@@ -2079 +2220 @@
                              mcevth[0].get_phi(), 
                              mcevth[0].get_core(),
-                             mcevth[0].get_coreX(),
-                             mcevth[0].get_coreY(),
+                             coreX,
+                             coreY,
                              impactD[0],
                              phiCT[0],
@@ -2116 +2257 @@
                              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 +2455 @@
               //
             
+
               for (int ict=0;ict<ct_Number;ict++){
                 Float_t ftime, ltime;
@@ -2327 +2469 @@
                                   mcevth[ict].get_phi(), 
                                   mcevth[ict].get_core(),
-                                  mcevth[ict].get_coreX(),
-                                  mcevth[ict].get_coreY(),
+                                  coreX,
+                                  coreY,
                                   impactD[ict],
                                   phiCT[ict],
@@ -2364 +2506 @@
                                   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 +2612 @@
           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 +2679 @@
   }
 
-  get_source_off(&sofftheta,&soffphi);
   if(!Trigger_Loop) icontrigger=0;
   time (&ltime);
@@ -2564 +2708 @@
                     telestheta,
                     telesphi,
-                    sofftheta,
-                    soffphi,
                     shthetamax,
                     shthetamin,
@@ -2605 +2747 @@
                     telestheta,
                     telesphi,
-                    sofftheta,
-                    soffphi,
                     shthetamax,
                     shthetamin,
@@ -2623 +2763 @@
   // Fill some missing values for MRawRunHeader
 
-  RunHeader->SetRunNumber(rnum);
+  RunHeader->SetRunNumber((UInt_t)rnum);
   RunHeader->SetNumEvents(nstoredevents);
 
@@ -2688 +2828 @@
     }
     
-    // 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 +3125 @@
       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 +3181 @@
 
   //------------------------------------------------------------
-  // Read Light Guides data
+  // Read Light Collection data
 
   // try to open the file
@@ -3055 +3200 @@
   // 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 +3218 @@
     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 +3719 @@
 // 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 +3732 @@
           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 +3807 @@
 int size_rotated(
     float *rotated,
-    float nsb[iMAXNUMPIX],
+    float *nsb,
     float rho)
 {
@@ -3737 +3908 @@
                   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 +3925 @@
   static float radius_mm;
 
+  Float_t t;
+
   // reset variables
 
@@ -3761 +3935 @@
   }
 
+  *itotnphe = 0;
   *incph = 0;
 
@@ -3768 +3943 @@
   random.SetSeed((Int_t)(get_seeds(0)*get_seeds(1)));
 
-  float spotsigma = get_sigma_xy_cm_spot();
+  float C1, C2, C3, rho;
+
 
   //- - - - - - - - - - - - - - - - - - - - - - - - - 
@@ -3779 +3955 @@
   
   // 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 +3993 @@
     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 +4046 @@
     
     // set pointer to the QE table of the relevant pixel
-          
+
     qept = (float **)QE[ict][ipixnum];
           
@@ -3870 +4067 @@
     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 +4111 @@
           
     nphe[ipixnum] += 1.0;
-          
+
+
     // store the new photoelectron
 
@@ -3929 +4139 @@
 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 +4170 @@
   char flag_new[4]; 
 
+
   // open input file 
   
@@ -3965 +4177 @@
   infile = fopen( iname, "r" );
 
-  // check if the satrfield input file exists
+  // check if the starfield input file exists
 
   if ( infile == NULL ) {
@@ -3975 +4187 @@
     return (0);
   }
-  
+
   // get signature, and check it
   
@@ -3982 +4194 @@
   }
 
-  // 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 +4263 @@
         
         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 +4308 @@
                               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 +4315 @@
                               &tmin_ns,
                               &tmax_ns,
-                              0
-                              ); // photons from starfield
+                              0,
+                              mirror_fraction); // photons from starfield
 
             if( k != 0 ){ // non-zero returnvalue means error
@@ -4124 +4351 @@
 
   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 +4363 @@
 //
 // $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 +4395 @@
 // 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.
 //