Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 2329)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 2334)
@@ -21,7 +21,7 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.57 $
+// $Revision: 1.58 $
 // $Author: blanch $ 
-// $Date: 2003-07-17 18:02:46 $
+// $Date: 2003-09-15 09:59:53 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -58,5 +58,6 @@
 
 #include "TArrayC.h"
-
+#include "TObjArray.h"
+                                                         
 #include "MTrigger.hxx"
 #include "MFadc.hxx"
@@ -75,4 +76,10 @@
 #include "MMcFadcHeader.hxx"
 #include "MGeomCamMagic.h"
+#include "MGeomCamMagic919.h"
+#include "MGeomCamMagicHG.h"
+#include "MGeomCamECO1000.h" 
+#include "MGeomCamECO1000HG.h" 
+#include "MGeomCamCT1.h"
+#include "MGeomCamCT1Daniel.h"
 #include "MGeomPix.h"
 
@@ -189,4 +196,10 @@
 static char ct_filename[256];  
 
+//@: Number of CT
+static int ct_Number;        
+
+//@: Camera geometries
+static int ct_Geometry;        
+
 //@: list of showers to be skipped
 static int *Skip;
@@ -245,7 +258,7 @@
 
 //@: Trigger conditions for a single trigger mode
-static float qThreshold[CAMERA_PIXELS];  
-static int Trigger_multiplicity = 4;
-static int Trigger_topology = 2;
+static float **qThreshold;  
+static int Trigger_multiplicity[100];
+static int Trigger_topology[100];
 
 //@: Upper and lower edges of the trigger loop
@@ -301,8 +314,8 @@
 
 //@: table of QE
-static float ***QE;
+static float ****QE;
 
 //@: number of datapoints for the QE curve
-static int pointsQE;
+static int pointsQE[100];
 
 //@: table of QE
@@ -398,5 +411,5 @@
   //@'
 
-  char inname[256];           //@< input file name
+  char inname_CT[100][256];   //@< array of names for each CT input file
   char starfieldname[256];    //@< starfield input file name
   char qe_filename[256];       //@< qe input file name
@@ -414,8 +427,11 @@
   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
-
-  int reflector_file_version;  //@< vrsion of he reflector file
-
-  FILE *inputfile;            //@< stream for the input file
+  int GeometryCamera[100];           //@< Identification of MGeomCam for each CT
+  int TriggerPixels[100];           //@< Numeber 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
+
   ofstream datafile;          //@< stream for the data file
 
@@ -426,4 +442,5 @@
 
   int inumphe;                //@< number of photoelectrons in an event from showers
+  int inumphe_CT[100];         //@< number of photoelectrons in an event from showers
   float inumphensb;             //@< number of photoelectrons in an event from nsb
 
@@ -442,4 +459,5 @@
   int ntshow=0;               //@< total number of showers
   int ncph=0;                 //@< partial number of photons in a given run
+  int ncph_system=0;          //@< partial number of photons in a given run
   int ntcph=0;                //@< total number of photons
 
@@ -454,10 +472,10 @@
   int nphe2NSB;               //@< From how many phe will we simulate NSB?
   float meanNSB;              //@< diffuse NSB mean value (phe per ns per central pixel)
-  float diffnsb_phepns[iMAXNUMPIX];  //@< diffuse NSB values for each pixel  
+  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[iMAXNUMPIX];                 //@< NSB values for each pixel
-  float nsb_phepns_rotated[iMAXNUMPIX];         //@< NSB values for each pixel after rotation
+  float nsb_phepns[100][iMAXNUMPIX];                 //@< NSB values for each pixel
+  float nsb_phepns_rotated[100][iMAXNUMPIX];         //@< NSB values for each pixel after rotation
 
   Float_t nsb_trigresp[TRIGGER_TIME_SLICES];    //@< array to write the trigger
@@ -479,10 +497,8 @@
   float zenfactor=1.0;  //  correction factor calculated from the extinction
 
-  int anaPixels;
-
   int nstoredevents = 0;
   int flagstoring = 0;
 
-  int ntrigger = 0;           //@< number of triggers in the whole file
+  int ntrigger[100];           //@< number of triggers in the whole file
   int btrigger = 0;           //@< trigger flag
   int ithrescount;            //@< counter for loop over threshold trigger
@@ -516,4 +532,5 @@
   Float_t soffphi = 0.0;
   UInt_t corsika = 5200 ;
+  Float_t maxpimpact = 0.0;
 
   // Star Field Rotation variables
@@ -522,8 +539,4 @@
   Float_t rho  , C3 , C2 , C1; 
 
-  MGeomCamMagic camgeom; // structure holding the camera definition
-
-  ct_NPixels=camgeom.GetNumPixels();
-
   //!@' @#### Definition of variables for |getopt()|.
   //@'
@@ -534,7 +547,13 @@
   //@'
 
+  qThreshold = new float *[100];
+  for(i=0;i<100;i++)
+    qThreshold[i] = new float [CAMERA_PIXELS]; 
+
   for(i=0;i<iMAXNUMPIX;i++){
-    nsb_phepns[i]=0;
-    nsb_phepns_rotated[i]=0;
+    for(int ict=0;ict<100;ict++){      
+      nsb_phepns[ict][i]=0;
+      nsb_phepns_rotated[ict][i]=0;
+    }
     for(j=0;j<iNUMWAVEBANDS;j++)
       nsbrate_phepns[i][j]=0.0;    //@< Starlight rates initialised at 0.0
@@ -597,9 +616,68 @@
   Data_From_STDIN = get_data_from_stdin();
 
+  // get number of telescopes and camera geometries
+
+  ct_Number=get_ct_number();
+  ct_Geometry=get_ct_geometry();
+
+  for(int i=0;i<ct_Number;i++){
+    ntrigger[i]=0;
+    if(ct_Geometry>=i*10){
+      GeometryCamera[i]=int(ct_Geometry/pow(10.0,i))%10;
+    }
+    else
+      GeometryCamera[i]=0;
+  }
+
+  // structure holding the camera definition
+
+  TObjArray camgeom;
+
+  for (int i=0; i<ct_Number;i++){
+    switch(GeometryCamera[i]){
+    case 1: camgeom[i]=new MGeomCamMagic;
+      TriggerPixels[i]=TRIGGER_PIXELS_1;
+      break;
+    case 2: camgeom[i]=new MGeomCamMagic919;
+      TriggerPixels[i]=TRIGGER_PIXELS_2;
+      break;
+    case 3: camgeom[i]=new MGeomCamMagicHG;
+      TriggerPixels[i]=TRIGGER_PIXELS_3;
+      break;
+    case 5: camgeom[i]=new MGeomCamECO1000;
+      TriggerPixels[i]=TRIGGER_PIXELS_5;
+      break;
+    case 6: camgeom[i]=new MGeomCamECO1000HG;
+      TriggerPixels[i]=TRIGGER_PIXELS_6;
+      break;
+    case 8: camgeom[i]=new MGeomCamCT1;
+      TriggerPixels[i]=TRIGGER_PIXELS_8;
+      break;
+    case 9: camgeom[i]=new MGeomCamCT1Daniel;
+      TriggerPixels[i]=TRIGGER_PIXELS_9;
+      break;
+    default: camgeom[i]=new MGeomCamMagic;
+      TriggerPixels[i]=TRIGGER_PIXELS_1;
+      break;
+    }
+  }
+  
+  ct_NPixels=0;
+
+  for(int ict=0;ict<ct_Number;ict++)
+    ct_NPixels=(((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels()>(UInt_t)ct_NPixels)?
+      (Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels():
+	ct_NPixels;
+
   // read WC and QE files
 
-  strcpy( qe_filename, get_qe_filename());
-
-  read_QE(qe_filename);
+  // FIX ME ! 
+  // Currently the PMT N of any camera with same average QE has the same QE.
+
+  QE = new float *** [ct_Number];
+  for(int ict=0;ict<ct_Number;ict++){
+    strcpy( qe_filename, get_qe_filename(ict));
+    read_QE(qe_filename,ict);
+  }
   read_WC();
 
@@ -627,4 +705,6 @@
 		       &FADC_resp_ampl_out, &FADC_resp_fwhm_out);
 
+  // FIXME --- tirgger properties may depend on the Camrea geometry
+
   get_Trigger_properties( &Trigger_gate_length, &Trigger_overlaping_time, &Trigger_response_ampl, &Trigger_response_fwhm);
 
@@ -635,15 +715,27 @@
     (Trigger_loop_umult-Trigger_loop_lmult+1)*
     (Trigger_loop_utop-Trigger_loop_ltop+1);
+
+  // Trigger loop operation is not implemented for Multi telscopes
+
+  if ( Trigger_Loop && ct_Number > 1 ){
+    cout<<"ERROR : camera::main : Trigger loop option is not";
+    cout<<" implemented for Multi Telescopes option. "<<icontrigger<<
+      " "<<ct_Number<<endl;
+    exit(1);
+  }
   
   if (!Trigger_Loop){
-    get_Trigger_Single (qThreshold, &Trigger_multiplicity, &Trigger_topology);
+    get_Trigger_Single (qThreshold, 
+			&Trigger_multiplicity[0], 
+			&Trigger_topology[0]);
     icontrigger=1;
   }
   else
-    get_threshold(qThreshold);
+    get_threshold(qThreshold[0]);
 
   // get filenames
 
-  strcpy( inname, get_input_filename() );
+  for(int ict=0;ict<ct_Number;ict++)
+    strcpy( inname_CT[ict], get_input_filename(ict) );
 
   strcpy( starfieldname, get_starfield_filename() );
@@ -678,5 +770,5 @@
       "%s:\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n\t%20s:\t%s\n",
       "Filenames",
-      "In", inname, 
+      "In", inname_CT[0], 
       "Stars", starfieldname,
       "NSB database (inner pixels)", nsbpathname,
@@ -702,7 +794,7 @@
    	"%s:\n\t%20s: %f\n\t%20s: %i\n\t%20s: %i\n",
     	"Single Trigger mode",
-    	"Threshold",qThreshold[0],
-    	"Multiplicity",Trigger_multiplicity,
-    	"Topology",Trigger_topology);
+    	"Threshold",qThreshold[0][0],
+    	"Multiplicity",Trigger_multiplicity[0],
+    	"Topology",Trigger_topology[0]);
   }   
   else{
@@ -711,6 +803,6 @@
     	"Single Trigger mode",
     	"Threshold","Different for each pixel",
-     	"Multiplicity",Trigger_multiplicity,
-    	"Topology",Trigger_topology);
+     	"Multiplicity",Trigger_multiplicity[0],
+    	"Topology",Trigger_topology[0]);
   }
   // log flags information
@@ -798,4 +890,5 @@
 
   Int_t Lev0, Lev1; 
+  Int_t Lev0MT[100], Lev1MT[100]; 
   
   fadcValues    =  new TArrayC(FADC_SLICES);
@@ -804,7 +897,4 @@
   // number of pixels for parameters
     
-  anaPixels = get_ana_pixels();
-  anaPixels = (anaPixels == -1) ? camgeom.GetNumPixels() : anaPixels;
-
   // Switched off writing TObject
 
@@ -818,23 +908,32 @@
   // initialise instance of Trigger and FADC classes
 
-  MTrigger  Trigger(Trigger_gate_length, 
+  MTrigger **Trigger_CT;
+  Trigger_CT = new MTrigger *[ct_Number];
+
+  for (int i=0; i<ct_Number;i++){
+    Trigger_CT[i] = new  MTrigger(TriggerPixels[i],
+		    ((MGeomCam*)(camgeom.UncheckedAt(i))),
+		    Trigger_gate_length, 
 		    Trigger_overlaping_time,
 		    Trigger_response_ampl, 
 		    Trigger_response_fwhm);  //@< A instance of the Class MTrigger 
-
+  }
+
+  for (int ict=0; ict<ct_Number;ict++){
   // Generate databse for trigger electronic noise
-  Trigger.SetElecNoise(Trigger_noise) ;
+    if (addElecNoise)    
+      Trigger_CT[ict]->SetElecNoise(Trigger_noise) ;
 
   // Set Right Discriminator threshold, taking into account trigger pixels
-  Trigger.CheckThreshold(&qThreshold[0]);
-
+    Trigger_CT[ict]->CheckThreshold(&qThreshold[ict][0],GeometryCamera[ict]);
+  }
   // Set flag in pixel 0 (not used for trigger) that indicates if secure pixel
   // is active: secureDiskThres*10000+riseDiskThres
 
   if(riseDiskThres>0.0)
-    qThreshold[0]=((UInt_t)secureDiskThres*100)*100+riseDiskThres;
+    qThreshold[0][0]=((UInt_t)secureDiskThres*100)*100+riseDiskThres;
   else
-    qThreshold[0]=0.0;
-
+    qThreshold[0][0]=0.0;
+  
   //  Initialise McTrig information class if we want to save trigger informtion
 
@@ -843,16 +942,22 @@
   MMcFadcHeader **HeaderFadc = NULL;
 
+  int numberBranches;
+
+  if (!Trigger_Loop)
+    numberBranches=ct_Number;
+  else
+    numberBranches=icontrigger;
 
   if (Write_McTrig){
 
-    McTrig = new MMcTrig * [icontrigger];
-  
-    for (i=0;i<icontrigger;i++) {
+    McTrig = new MMcTrig * [numberBranches];
+  
+    for (i=0;i<numberBranches;i++) {
       McTrig[i] = new MMcTrig();
     }
 
-    HeaderTrig = new MMcTrigHeader * [icontrigger];
-  
-    for (i=0;i<icontrigger;i++) {
+    HeaderTrig = new MMcTrigHeader * [numberBranches];
+  
+    for (i=0;i<numberBranches;i++) {
       HeaderTrig[i] = new MMcTrigHeader();
     }
@@ -861,13 +966,20 @@
   if (Write_McFADC){
 
-    HeaderFadc = new MMcFadcHeader * [icontrigger];
-  
-    for (i=0;i<icontrigger;i++) {
+    HeaderFadc = new MMcFadcHeader * [numberBranches];
+  
+    for (i=0;i<numberBranches;i++) {
       HeaderFadc[i] = new MMcFadcHeader();
     }
   }
 
-  MFadc fadc(FADC_response_ampl,FADC_response_fwhm,
-	     FADC_resp_ampl_out,FADC_resp_fwhm_out) ; //@< A instance of the Class MFadc
+  MFadc **Fadc_CT;
+  Fadc_CT = new MFadc *[ct_Number];
+
+  for (int i=0; i<ct_Number;i++)
+    Fadc_CT[i] = 
+      new  MFadc(((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels(),
+		 FADC_response_ampl,FADC_response_fwhm,
+		 FADC_resp_ampl_out,FADC_resp_fwhm_out) ; //@< A instance of the Class MFadc
+
 
   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -883,24 +995,33 @@
   /////////////////////////////////////////////////////////////////////
 
-  Float_t input_pedestals[CAMERA_PIXELS];
-
-  for(i=0;i<CAMERA_PIXELS;i++)
+  Float_t input_pedestals[ct_NPixels];
+
+  for(i=0;i<ct_NPixels;i++)
     input_pedestals[i]=get_FADC_pedestal();
-
-  fadc.SetPedestals(input_pedestals);
-
-  // Generate databse for the Fadc electronic noise
-  
-  fadc.SetElecNoise(FADC_noise);
-
+  for (int ict=0; ict<ct_Number;ict++){
+
+    Fadc_CT[ict]->SetPedestals(input_pedestals);
+    
+    // Generate databse for the Fadc electronic noise
+    if (addElecNoise)        
+      Fadc_CT[ict]->SetElecNoise(FADC_noise);
+  }
   // Prepare the raw data output
 
        // Header Tree
-
+  MGeomCam **camdummy = new MGeomCam * [ct_Number];
+  for (int ict=0; ict<ct_Number;ict++){
+    camdummy[ict]= ((MGeomCam*)(camgeom.UncheckedAt(ict)));
+  }
   MRawRunHeader *RunHeader= new MRawRunHeader();
   MMcRunHeader  *McRunHeader = new MMcRunHeader();
-  MMcConfigRunHeader  *McConfigRunHeader = new MMcConfigRunHeader();
   MMcCorsikaRunHeader  *McCorsikaRunHeader = new MMcCorsikaRunHeader();
 
+  MMcConfigRunHeader  **McConfigRunHeader = NULL;
+  McConfigRunHeader = new MMcConfigRunHeader * [numberBranches];
+  for (i=0;i<numberBranches;i++) {
+    McConfigRunHeader[i] = new MMcConfigRunHeader();
+  }
+
        // Header branch
 
@@ -908,7 +1029,7 @@
 
   if (Write_RawEvt) {
-    EvtHeader = new MRawEvtHeader * [icontrigger]; 
-
-    for (i=0;i<icontrigger;i++) {
+    EvtHeader = new MRawEvtHeader * [numberBranches]; 
+
+    for (i=0;i<numberBranches;i++) {
       EvtHeader[i] = new MRawEvtHeader();
     }
@@ -920,7 +1041,7 @@
 
   if (Write_RawEvt) {
-    EvtData = new MRawEvtData * [icontrigger]; 
-
-    for (i=0;i<icontrigger;i++) {
+    EvtData = new MRawEvtData * [numberBranches]; 
+
+    for (i=0;i<numberBranches;i++) {
       EvtData[i] = new MRawEvtData();
       EvtData[i]->Init(RunHeader);     //  We need thr RunHeader to read
@@ -950,17 +1071,29 @@
   		    &McRunHeader);
 
-  HeaderTree.Branch("MMcConfigRunHeader","MMcConfigRunHeader",
-  		    &McConfigRunHeader);
-
   HeaderTree.Branch("MMcCorsikaRunHeader","MMcCorsikaRunHeader",
   		    &McCorsikaRunHeader);
 
-  if(!Trigger_Loop && Write_McTrig){
+  if(!Trigger_Loop && Write_McTrig && ct_Number==1){
     
     HeaderTree.Branch("MMcTrigHeader","MMcTrigHeader", 
 		      &HeaderTrig[0]);
   }
-  if (Trigger_Loop && Write_McTrig){
-    for(char branchname[20],i=0;i<icontrigger;i++){
+  if (ct_Number==1)
+    HeaderTree.Branch("MGeomCam","MGeomCam",
+		      &camdummy[0]);
+  else{
+    char branchname[20];
+    for (int ict=0; ict<ct_Number;ict++){
+      sprintf(help,"%i",ict+1);
+      strcpy (branchname, "MGeomCam;");
+      strcat (branchname, & help[0]);
+      strcat (branchname, ".");
+      HeaderTree.Branch(branchname,"MGeomCam", 
+			&camdummy[ict]);
+    }
+  }
+
+  if ((Trigger_Loop || ct_Number>1) && Write_McTrig){
+    for(char branchname[20],i=0;i<numberBranches;i++){
       
       sprintf(help,"%i",i+1);
@@ -973,11 +1106,11 @@
   }  
 
-  if(!Trigger_Loop && Write_McFADC){
+  if(!Trigger_Loop && Write_McFADC && ct_Number==1){
     
     HeaderTree.Branch("MMcFadcHeader","MMcFadcHeader", 
 		      &HeaderFadc[0]);
   }
-  if (Trigger_Loop && Write_McFADC){
-    for(char branchname[20],i=0;i<icontrigger;i++){
+  if ((Trigger_Loop || ct_Number>1) && Write_McFADC){
+    for(char branchname[20],i=0;i<numberBranches;i++){
       
       sprintf(help,"%i",i+1);
@@ -999,14 +1132,13 @@
   RunHeader->SetNumSamples(0,FADC_SLICES);
   RunHeader->SetNumCrates(1);
-  RunHeader->SetNumPixInCrate(CAMERA_PIXELS);
-
+  RunHeader->SetNumPixInCrate(ct_NPixels);
 
   //  Fill branches for MMcTrigHeader
 
-  if(!Trigger_Loop && Write_McTrig){
-
-    HeaderTrig[0]->SetTopology((Short_t) Trigger_topology);
-    HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity);
-    HeaderTrig[0]->SetThreshold(qThreshold);
+  if(!Trigger_Loop && Write_McTrig && ct_Number==1){
+
+    HeaderTrig[0]->SetTopology((Short_t) Trigger_topology[0]);
+    HeaderTrig[0]->SetMultiplicity((Short_t) Trigger_multiplicity[0]);
+    HeaderTrig[0]->SetThreshold(qThreshold[0]);
     HeaderTrig[0]->SetAmplitud(Trigger_response_ampl);
     HeaderTrig[0]->SetFwhm(Trigger_response_fwhm);
@@ -1015,4 +1147,16 @@
     HeaderTrig[0]->SetElecNoise(Trigger_noise);
   }
+  if(!Trigger_Loop && Write_McTrig && ct_Number>1){
+    for(int i=0;i<ct_Number;i++){
+      HeaderTrig[i]->SetTopology((Short_t) Trigger_topology[i]);
+      HeaderTrig[i]->SetMultiplicity((Short_t) Trigger_multiplicity[i]);
+      HeaderTrig[i]->SetThreshold(qThreshold[i]);
+      HeaderTrig[i]->SetAmplitud(Trigger_response_ampl);
+      HeaderTrig[i]->SetFwhm(Trigger_response_fwhm);
+      HeaderTrig[i]->SetOverlap(Trigger_overlaping_time);
+      HeaderTrig[i]->SetGate(Trigger_gate_length);
+      HeaderTrig[i]->SetElecNoise(Trigger_noise);
+    }    
+  }  
   if(Trigger_Loop && Write_McTrig){
 
@@ -1023,8 +1167,8 @@
 	  HeaderTrig[iconcount]->SetTopology((Short_t) isorttopo[itopocount+Trigger_loop_ltop]);
 	  HeaderTrig[iconcount]->SetMultiplicity((Short_t) imulticount+Trigger_loop_lmult);
-	  for(int i=0;i<CAMERA_PIXELS;i++){
+	  for(int i=0;i<ct_NPixels;i++){
 	    fpixelthres[i]=
-	      ((Float_t)(fthrescount)>=qThreshold[i])?
-	      (Float_t)(fthrescount):qThreshold[i];
+	      ((Float_t)(fthrescount)>=qThreshold[0][i])?
+	      (Float_t)(fthrescount):qThreshold[0][i];
 	      }
 	  HeaderTrig[iconcount]->SetThreshold( fpixelthres);
@@ -1039,21 +1183,32 @@
     }
   }
-
+  
   //  Fill branches for MMcFadcHeader
 
-  for(i=0;i<CAMERA_PIXELS;i++){
+  for(i=0;i<ct_NPixels;i++){
       fadc_elecnoise[i]=FADC_noise;
   }
 
-  fadc.GetPedestals(&fadc_pedestals[0]);
-
-  if(!Trigger_Loop && Write_McFADC){
+  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);
     HeaderFadc[0]->SetFwhm(FADC_response_fwhm);
-    HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],anaPixels);
-    HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0],anaPixels);
-  }
+    HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],
+			       ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
+    HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0],
+				((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
+  }
+  if(!Trigger_Loop && Write_McFADC && ct_Number>1){
+    for(int i=0;i<ct_Number;i++){
+      HeaderFadc[i]->SetShape(0.0);
+      HeaderFadc[i]->SetAmplitud(FADC_response_ampl);
+      HeaderFadc[i]->SetFwhm(FADC_response_fwhm);
+      HeaderFadc[i]->SetPedestal(&fadc_pedestals[0],((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
+      HeaderFadc[i]->SetElecNoise(&fadc_elecnoise[0],((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
+    }
+  }  
   if(Trigger_Loop && Write_McFADC){
     int iconcount;
@@ -1064,6 +1219,6 @@
 	  HeaderFadc[iconcount]->SetAmplitud(Trigger_response_ampl);
 	  HeaderFadc[iconcount]->SetFwhm(Trigger_response_fwhm);
-	  HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0],anaPixels);
-	  HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],anaPixels);
+	  HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0], ct_NPixels);
+	  HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0], ct_NPixels);
 	  iconcount++;
 	}
@@ -1072,7 +1227,4 @@
   }
 
-  //  Fill the Header Tree with the current leaves of each branch
-  // HeaderTree.Fill() ;
-	    
 
   //      create a Tree for the Event data stream 
@@ -1085,6 +1237,9 @@
   }
 
-  if(!Trigger_Loop){
+  if(!Trigger_Loop && ct_Number==1){
     
+    HeaderTree.Branch("MMcConfigRunHeader","MMcConfigRunHeader",
+		      &McConfigRunHeader[0]);
+
     if (Write_RawEvt){
       EvtTree.Branch("MRawEvtHeader","MRawEvtHeader", 
@@ -1099,6 +1254,18 @@
   }
   else{
+
+    if(ct_Number>1)
+      for(char branchname[10],i=0;i<numberBranches;i++){
+      
+	sprintf(help,"%i",i+1);
+	strcpy (branchname, "MMcConfigRunHeader;");
+	strcat (branchname, & help[0]);
+	strcat (branchname, ".");
+	HeaderTree.Branch(branchname,"MMcConfigRunHeader",
+			  &McConfigRunHeader[i]);
+      }
+    
     if (Write_McTrig){
-      for(char branchname[10],i=0;i<icontrigger;i++){
+      for(char branchname[10],i=0;i<numberBranches;i++){
       
 	sprintf(help,"%i",i+1);
@@ -1112,6 +1279,6 @@
   }  
 
-  if (Trigger_Loop && Write_RawEvt){
-    for(char branchname[15],i=0;i<icontrigger;i++){
+  if ((Trigger_Loop || ct_Number>1) && Write_RawEvt){
+    for(char branchname[15],i=0;i<numberBranches;i++){
       
       sprintf(help,"%i",i+1);
@@ -1122,5 +1289,5 @@
 		     &EvtHeader[i]);
     }
-    for(char branchname[15],i=0;i<icontrigger;i++){
+    for(char branchname[15],i=0;i<numberBranches;i++){
       
       sprintf(help,"%i",i+1);
@@ -1133,4 +1300,5 @@
   }
 
+
   TApplication theAppTrigger("App", &argc, argv);
 
@@ -1171,7 +1339,9 @@
     //
     
+    // FIXME --- star NSB different for each camera?
     k = produce_nsbrates( starfieldname,
-			  &camgeom,
-			  nsbrate_phepns);
+			  ((MGeomCam*)(camgeom.UncheckedAt(0))),
+			  nsbrate_phepns,
+			  0);
 
     if (k != 0){
@@ -1180,25 +1350,56 @@
     }
   
-    // calculate diffuse rate correcting for the pixel size
-
-    const double size_inner = camgeom[0].GetR();
-      
-    for(UInt_t ui=0; ui<camgeom.GetNumPixels(); ui++){
-      const double actual_size = camgeom[ui].GetR();
-      diffnsb_phepns[ui] = meanNSB * 
-	actual_size*actual_size/(size_inner*size_inner);
+    // calculate diffuse rate correcting for the pixel size and telescope
+
+    float factorNSB_ct;
+    for(int ict=0;ict<ct_Number;ict++){      
+      
+      // First we set the factor due to the mirror size respect to the MAGIC mirror.
+      
+      switch(GeometryCamera[ict]){
+      case 1:
+      case 2:
+      case 3:
+      case 4:
+	factorNSB_ct=1.0;
+	break;
+      case 5:
+      case 6:
+      case 7:
+	factorNSB_ct=1000.0/239.0;
+	break;
+      case 8:
+      case 9:
+      default:
+	factorNSB_ct=1.0;
+	break;	
+      }
+
+      factorNSB_ct=factorNSB_ct/
+	((*((MGeomCam*)(camgeom.UncheckedAt(ict)))).GetCameraDist()*1000*
+	 (*((MGeomCam*)(camgeom.UncheckedAt(ict)))).GetCameraDist()*1000*
+	 PIXEL_SIZE*PIXEL_SIZE);
+
+      for(UInt_t ui=0;
+	  ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels(); ui++){
+	diffnsb_phepns[ict][ui] = (Int_t(meanNSB*factorNSB_ct*100+0.5))/100.0 * 
+	  (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD()*
+	  (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD();
+      }
     }
 
     // calculate nsb rate including diffuse and starlight
     // we also include the extinction effect
-    for(j=0;j<iNUMWAVEBANDS;j++){
+    for(int ict=0;ict<ct_Number;ict++){      
+      for(j=0;j<iNUMWAVEBANDS;j++){
 	// calculate the effect of the atmospheric extinction 
-    
+	
 	zenfactor = pow(10., -0.4 * ext[j] );
 	
-	for(UInt_t ui=0; ui<camgeom.GetNumPixels();ui++){
-	    nsb_phepns[ui]+=diffnsb_phepns[ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
-	    nsb_phepns_rotated[ui]=nsb_phepns[ui];
+	for(UInt_t ui=0; ui<((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels();ui++){
+	  nsb_phepns[ict][ui]+=diffnsb_phepns[ict][ui]/iNUMWAVEBANDS + zenfactor *nsbrate_phepns[ui][j];
+	  nsb_phepns_rotated[ict][ui]=nsb_phepns[ict][ui];
 	}
+      }
     }
   }
@@ -1212,21 +1413,27 @@
   if ( Data_From_STDIN ) {
 
-    inputfile = stdin;
+    inputfile[0] = stdin;
 
   }
   else{
 
-    log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname );
-    inputfile = fopen( inname, "r" );
-    if ( inputfile == NULL ) 
-      error( SIGNATURE, "Cannot open input file: %s\n", inname );
-
+    for(int ict=0;ict<ct_Number;ict++){
+      
+      log( SIGNATURE, "Opening input \"rfl\" file %s\n", inname_CT[ict] );
+      inputfile[ict] = fopen( inname_CT[ict], "r" );
+
+      if ( inputfile[ict] == NULL ) 
+	error( SIGNATURE, "Cannot open input file: %s\n", inname_CT[ict] );
+    }
+    
   }
   
   // get signature, and check it
 
-  if((reflector_file_version=check_reflector_file( inputfile ))==FALSE){
-    exit(1);
-  }
+  for(int ict=0;ict<ct_Number;ict++){
+    if((reflector_file_version=check_reflector_file( inputfile[ict] ))==FALSE){
+      exit(1);
+    }
+  }  
 
   // open data file
@@ -1245,5 +1452,5 @@
   // allocate space for PMTs numbers of pixels
   
-  fnpix = new float [ camgeom.GetNumPixels() ];
+  fnpix = new float [ct_NPixels];
 
      
@@ -1252,6 +1459,8 @@
 
   // get flag
-    
-  fread( flag, SIZE_OF_FLAGS, 1, inputfile );
+
+  for(int ict=0;ict<ct_Number;ict++)   
+    fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
+  
 
   // loop over the file
@@ -1259,6 +1468,8 @@
   still_in_loop = TRUE;
 
+  // FIXME --- check if not eof for all input reflector files
+
   while (
-         ((! Data_From_STDIN) && ( !feof(inputfile) ))
+         ((! Data_From_STDIN) && ( !feof(inputfile[0]) ))
          ||
          (Data_From_STDIN && still_in_loop)
@@ -1275,41 +1486,47 @@
     else { // found start of run
 
-      fread( flag_new, 4, 1, inputfile );
-
-      if(!isA( flag_new, FLAG_START_OF_HEADER)){
-
-	//  We break the main loop
-	cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
-	cout<<"         We break the main loop"<<endl;
-	break;      
-      }
-
-      //Float_t flag_temp[SIZE_OF_MCRUNHEADER];
-      
-      //fread( flag_temp, (SIZE_OF_MCRUNHEADER)*sizeof(float), 1, inputfile );
-      fread((char*)&mcrunh,(SIZE_OF_MCRUNHEADER)*sizeof(float), 1, inputfile );
-
-      nshow=0;
-
-      fread( flag, SIZE_OF_FLAGS, 1, inputfile );
-      while( isA( flag, FLAG_START_OF_EVENT   )){ // while there is a next event
-	fread( flag_new, 4, 1, inputfile );
-
-	if(!isA( flag_new, FLAG_EVENT_HEADER)){
-
-	//  We break while events loop
-	  cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<endl;
-	  cout<<"         We break the while events loop"<<endl;
+      for(int ict=0;ict<ct_Number;ict++) {
+	
+	fread( flag_new, 4, 1, inputfile[ict] );
+
+	if(!isA( flag_new, FLAG_START_OF_HEADER)){
+	  
+	  //  We break the main loop
+	  cout<<"Warning: Expected start of run header flag, but found:"<<flag_new<<endl;
+	  cout<<"         We break the main loop"<<endl;
 	  break;      
 	}
+
+	fread((char*)&mcrunh,(SIZE_OF_MCRUNHEADER)*sizeof(float), 1, inputfile[ict] );
+      }
+      
+      nshow=0;
+
+      for(int ict=0;ict<ct_Number;ict++){ 
+	fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
+      }
+      while( isA( flag, FLAG_START_OF_EVENT   )){ // while there is a next event
+	for(int ict=0;ict<ct_Number;ict++) {
+	  fread( flag_new, 4, 1, inputfile[ict] );
+	  if(!isA( flag_new, FLAG_EVENT_HEADER)){
+
+	    //  We break while events loop
+	    cout<<"Warning: Expected start of event header flag, but found:"<<flag_new<<" "<<ict<<endl;
+	    cout<<"         We break the while events loop"<<endl;
+	    break;      
+	  }
+	}
+	
 
 	//
 	// Clear Trigger and Fadc
 	//
-      	Trigger.Reset() ; 
-	Trigger.ClearFirst();
-	Trigger.ClearZero();
-	fadc.Reset() ; 
-
+	for(int ict=0;ict<ct_Number;ict++) {
+	  Trigger_CT[ict]->Reset() ; 
+	  Trigger_CT[ict]->ClearFirst();
+	  Trigger_CT[ict]->ClearZero();
+	  Fadc_CT[ict]->Reset() ; 
+	}
+	
 	++nshow;
 	if((nshow+ntshow)%100 == 1)
@@ -1318,8 +1535,10 @@
 	// get MCEventHeader
 	if (reflector_file_version<6){
-	  fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile );
+	  for(int ict=0;ict<ct_Number;ict++)
+	    fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile[ict] );
 	}
 	else{
-	  fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile );
+	  for(int ict=0;ict<ct_Number;ict++)
+	    fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile[ict] );
 	}
 
@@ -1461,5 +1680,5 @@
 	}
 
-	// Read first and last time and put inumphe to 0
+	// Read first and last time and put inumphe_CT[0] to 0
 
 	if (reflector_file_version<6)
@@ -1468,26 +1687,34 @@
 	  mcevth_2.get_times(&arrtmin_ns,&arrtmax_ns);
 
+	ncph_system=0;
 	inumphe=0;
 
-	// read the photons and produce the photoelectrons
-
-	k = produce_phes( inputfile,
-			  &camgeom,
-			  WAVEBANDBOUND1,
-			  WAVEBANDBOUND6,
-			  &Trigger,   // will be changed by the function!
-			  &fadc,      // will be changed by the function!
-			  &inumphe, // important for later: the size of photoe[]
-			  fnpix,    // will be changed by the function!
-			  &ncph,    // will be changed by the function!
-			  &arrtmin_ns, // will be changed by the function!
-			  &arrtmax_ns); // will be changed by the function!
-
-
-	if( k != 0 ){ // non-zero returnvalue means error
-	  cout << "Exiting.\n";
-	  exit(1);
+	for(int ict=0;ict<ct_Number;ict++){
+	  
+	  inumphe_CT[ict]=0;
+
+	  // read the photons and produce the photoelectrons
+	  
+	  k = produce_phes( inputfile[ict],
+			    ((MGeomCam*)(camgeom.UncheckedAt(ict))),
+			    WAVEBANDBOUND1,
+			    WAVEBANDBOUND6,
+			    Trigger_CT[ict],   // will be changed by the function!
+			    Fadc_CT[ict],      // will be changed by the function!
+			    &inumphe_CT[ict], // important for later: the size of photoe[]
+			    fnpix,    // will be changed by the function!
+			    &ncph,    // will be changed by the function!
+			    &arrtmin_ns, // will be changed by the function!
+			    &arrtmax_ns, // will be changed by the function!
+			    ict);
+
+	  inumphe=(inumphe<inumphe_CT[ict])?inumphe_CT[ict]:inumphe;
+	  
+	  if( k != 0 ){ // non-zero returnvalue means error
+	    cout << "Exiting.\n";
+	    exit(1);
+	  }
 	}
-
+	
 	// NSB simulation
 
@@ -1516,43 +1743,49 @@
 
 	    // Rottion of the NSB
-		
-	    k = size_rotated(
-			     &nsb_phepns_rotated[0],
-			     nsb_phepns,
-			     rho);
+	    // FIXME --- We should rotate for all cameras. Is it always the same rho? 	
+	    for(int ict=0;ict<ct_Number;ict++){      
+	      k = size_rotated(
+			       &nsb_phepns_rotated[ict][0],
+			       nsb_phepns[ict],
+			       rho);
+	    }
+	  }
+
+	  //  Fill trigger and fadc response in the trigger class from the database
+	  for(int ict=0;ict<ct_Number;ict++){
 	    
+	    for(UInt_t ui=0;ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++)
+	      if(nsb_phepns_rotated[ict][ui]>0.0){
+		if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD()>(*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD()){
+		  k=lons_outer.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
+					   & nsb_trigresp[0],
+					   & nsb_fadcresp[0]);
+	      }
+		else{
+		  k=lons.GetResponse(nsb_phepns_rotated[ict][ui],0.01,
+				     & nsb_trigresp[0],& nsb_fadcresp[0]);
+		}
+		if(k==0){
+		  cout << "Exiting.\n";
+		  exit(1);
+		}
+		Trigger_CT[ict]->AddNSB(ui,nsb_trigresp);
+		Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
+	      }
 	  }
-
-	  //  Fill trigger and fadc response in the trigger class from the database
-	  for(UInt_t ui=0;ui<camgeom.GetNumPixels();ui++)
-	    if(nsb_phepns_rotated[ui]>0.0){
-	      if(camgeom[ui].GetR()>camgeom[0].GetR()){
-		k=lons_outer.GetResponse(nsb_phepns_rotated[ui],0.1,
-				   & nsb_trigresp[0],& nsb_fadcresp[0]);
-	      }
-	      else{
-		k=lons.GetResponse(nsb_phepns_rotated[ui],0.1,
-				   & nsb_trigresp[0],& nsb_fadcresp[0]);
-	      }
-	      if(k==0){
-		cout << "Exiting.\n";
-		exit(1);
-	      }
- 	      Trigger.AddNSB(ui,nsb_trigresp);
-	      fadc.AddSignal(ui,nsb_fadcresp);
-	    }
-	}// end if(simulateNSB && nphe2NSB<=inumphe) ...
+	
+	}// end if(simulateNSB && nphe2NSB<=inumphe_CT[0]) ...
   
 	inumphensb=0;
-	for (UInt_t ui=0;ui<camgeom.GetNumPixels();ui++){
-	  inumphensb+=nsb_phepns[ui]*TOTAL_TRIGGER_TIME;
+	for (UInt_t ui=0;ui<
+	       ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels();ui++){
+	  inumphensb+=nsb_phepns[0][ui]*TOTAL_TRIGGER_TIME;
 	}
-
-	ntcph+=ncph;
+	ntcph+=ncph_system;
 	if ((nshow+ntshow)%100 == 1){
 	  log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
 	      ncph, ntcph);
 
-	  cout << "Total number of phes: " << inumphe<<" + ";
+	  cout << "Total number of phes in CT 0: " << inumphe_CT[0]<<" + ";
 	  cout<<inumphensb<<endl;
 	}
@@ -1600,15 +1833,20 @@
 	//   
 
-	if (addElecNoise && nphe2NSB<=inumphe){	
-	  Trigger.ElecNoise(Trigger_noise) ;
+	for(int ict=0;ict<ct_Number;ict++) {
+	  if (addElecNoise && nphe2NSB<=inumphe){	
+
+	    Trigger_CT[ict]->ElecNoise(Trigger_noise) ;
 	    
-	  fadc.ElecNoise(FADC_noise) ;
+	    Fadc_CT[ict]->ElecNoise(FADC_noise) ;
+	  }
 	}
-
+	
 	//  now a shift in the fadc signal due to the pedestlas is 
 	//  intorduced
 	//  This is done inside the class MFadc by the method Pedestals
-	fadc.Pedestals();
-
+	for(int ict=0;ict<ct_Number;ict++) {
+	  Fadc_CT[ict]->Pedestals();
+	}
+	
 	//   We study several trigger conditons
 	if(Trigger_Loop){
@@ -1617,22 +1855,25 @@
 	  btrigger=0;
 	  flagstoring = 0;
+
 	  //  Loop over trigger threshold
 	  int iconcount;
 	  for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
-	    for (i=0;i<CAMERA_PIXELS;i++){
+	    for (i=0;i<ct_NPixels;i++){
 	      fpixelthres[i]=
-		((Float_t)(fthrescount)>=qThreshold[i])?
-		(Float_t)(fthrescount):qThreshold[i];
+		((Float_t)(fthrescount)>=qThreshold[0][i])?
+		(Float_t)(fthrescount):qThreshold[0][i];
 
 	      // Rise the discrimnator threshold to avoid huge rates
 
 	      if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
-		for(int ii=0;ii<CAMERA_PIXELS;ii++)
-		  if( nsb_phepns_rotated[ii]>riseDiskThres)
+		for(int ii=0;ii<ct_NPixels;ii++)
+		  if( nsb_phepns_rotated[0][ii]>riseDiskThres)
 		    fpixelthres[ii]=secureDiskThres;
 	    }
-	    Trigger.SetThreshold(fpixelthres);
-
-	    Trigger.Diskriminate();
+	    Trigger_CT[0]->SetThreshold(fpixelthres);
+	    
+	    Trigger_CT[0]->Diskriminate();
+
+	    
 	    //
 	    //   look if in all the signals in the trigger signal branch
@@ -1645,12 +1886,13 @@
 	    
 	    //  Set trigger flags to zero
+	    Lev0=0;
 	    Lev1=0;
 
 	    //  loop over multiplicity of trigger configuration
 	    for (imulticount=Trigger_loop_lmult;imulticount<=Trigger_loop_umult;imulticount++){
-	      Trigger.SetMultiplicity(imulticount);
-	      Trigger.ClearZero();
-
-	      Lev0=(Short_t) Trigger.ZeroLevel();
+	      Trigger_CT[0]->SetMultiplicity(imulticount);
+	      Trigger_CT[0]->ClearZero();
+
+	      Lev0=(Short_t) Trigger_CT[0]->ZeroLevel();
 	      if (Lev0>0 || Write_All_Images || btrigger){
 		
@@ -1664,12 +1906,12 @@
 		  // if there are 3 or more N pixels. 
 		  if(imulticount<3)
-		    Trigger.SetTopology(1);
+		    Trigger_CT[0]->SetTopology(1);
 		  else
 		    {
 		      // We should be careful that topologies are sort from 
 		      // the less to the more restrictive one. 
-		      Trigger.SetTopology(isorttopo[itopocount]);
+		      Trigger_CT[0]->SetTopology(isorttopo[itopocount]);
 		    }
-		  Trigger.ClearFirst();
+		  Trigger_CT[0]->ClearFirst();
 		  
 		  //
@@ -1677,5 +1919,5 @@
 		  //
 		  if(Lev0!=0)
-		      Lev1=Trigger.FirstLevel();
+		      Lev1=Trigger_CT[0]->FirstLevel();
 		  if(Lev1>0) {
 		    btrigger= 1;
@@ -1694,6 +1936,6 @@
 		      if (Write_McTrig){
 			  McTrig[iconcount]->SetFirstLevel ((ii+1)*Lev0);
-			  McTrig[iconcount]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
-			  Trigger.GetMapDiskriminator(trigger_map);
+			  McTrig[iconcount]->SetTime(Trigger_CT[0]->GetFirstLevelTime(ii),ii+1);
+			  Trigger_CT[0]->GetMapDiskriminator(trigger_map);
 			  McTrig[iconcount]->SetMapPixels(trigger_map,ii);
 		      }
@@ -1702,5 +1944,5 @@
 		    //
 
-		    fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
+		    Fadc_CT[0]->TriggeredFadc(Trigger_CT[0]->GetFirstLevelTime(ii));
 	    
 		    if( Write_RawEvt ){
@@ -1712,8 +1954,10 @@
 		      //   fill pixel information
 		      if (Lev1){
-			for(UInt_t i=0;i<camgeom.GetNumPixels();i++){
+			for(UInt_t i=0;
+			    i<((MGeomCam*)(camgeom.UncheckedAt(0)))
+			      ->GetNumPixels();i++){
 			  for (j=0;j<FADC_SLICES;j++){
-			    fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
-			    fadcValuesLow->AddAt(fadc.GetFadcLowGainSignal(i,j),j);
+			    fadcValues->AddAt(Fadc_CT[0]->GetFadcSignal(i,j),j);
+			    fadcValuesLow->AddAt(Fadc_CT[0]->GetFadcLowGainSignal(i,j),j);
 			  }
 			  EvtData[iconcount]->AddPixel(i,fadcValues,0);
@@ -1776,6 +2020,6 @@
 			     (UInt_t)(mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot()), 
 			     (UInt_t) ncph, 
-			     (UInt_t) inumphe,
-			     (UInt_t) inumphensb+inumphe,
+			     (UInt_t) inumphe_CT[0],
+			     (UInt_t) inumphensb+inumphe_CT[0],
 			     -1.0,
 			     -1.0,
@@ -1813,6 +2057,6 @@
 			     (UInt_t)(mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot()), 
 			     (UInt_t) ncph, 
-			     (UInt_t) inumphe,
-			     (UInt_t) inumphensb+inumphe,
+			     (UInt_t) inumphe_CT[0],
+			     (UInt_t) inumphensb+inumphe_CT[0],
 			     mcevth_2.get_ElecFraction(),
 			     mcevth_2.get_MuonFraction(),
@@ -1825,10 +2069,10 @@
 	    //  Clear the branches
 	    if(Write_McTrig){
-	      for(i=0;i<icontrigger;i++){
+	      for(i=0;i<numberBranches;i++){
 		McTrig[i]->Clear() ;
 	      }
 	    }
 	    if( Write_RawEvt ){
-	      for(i=0;i<icontrigger;i++){
+	      for(i=0;i<numberBranches;i++){
 		EvtHeader[i]->Clear() ;
 		EvtData[i]->DeletePixels();
@@ -1842,91 +2086,142 @@
 	else {
 
-	  //  Setting trigger conditions
-	  Trigger.SetMultiplicity(Trigger_multiplicity);
-	  Trigger.SetTopology(Trigger_topology);
-	  for (int i=0;i<CAMERA_PIXELS;i++)
-	    fpixelthres[i]=qThreshold[i];
-
-	  // Rise the discrimnator threshold to avoid huge rates
-	  if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
-	    for(int ii=0;ii<CAMERA_PIXELS;ii++){
-	      if( nsb_phepns_rotated[ii]>riseDiskThres)
-		fpixelthres[ii]=secureDiskThres;
+	  // Set to zero the flag to know if some conditon has triggered
+	  btrigger=0;
+	  flagstoring = 0;
+	  
+	  for(int ict=0;ict<ct_Number;ict++){
+
+	    //  Setting trigger conditions
+	    Trigger_CT[ict]->SetMultiplicity(Trigger_multiplicity[ict]);
+	    Trigger_CT[ict]->SetTopology(Trigger_topology[ict]);
+	    for (int i=0;i<ct_NPixels;i++)
+	      fpixelthres[i]=qThreshold[ict][i];
+
+	    // Rise the discrimnator threshold to avoid huge rates
+	    if(riseDiskThres>0.0 && simulateNSB && nphe2NSB<=inumphe)
+	      for(int ii=0;ii<ct_NPixels;ii++){
+		if( nsb_phepns_rotated[ict][ii]>riseDiskThres)
+		  fpixelthres[ii]=secureDiskThres;
+	      }
+	  
+	    Trigger_CT[ict]->SetThreshold(fpixelthres);
+				 
+	    Trigger_CT[ict]->Diskriminate() ; 
+
+	    //
+	    //   look if in all the signals in the trigger signal branch
+	    //   is a possible Trigger. Therefore we habe to diskriminate all
+	    //   the simulated analog signals (Method Diskriminate in class
+	    //   MTrigger). We look simultanously for the moments at which
+	    //   there are more than TRIGGER_MULTI pixels above the 
+	    //   CHANNEL_THRESHOLD. 
+	    //
+	    
+	    Lev0MT[ict] = (Short_t) Trigger_CT[ict]->ZeroLevel() ; 
+	  
+	    Lev1MT[ict] = 0 ; 
+	  
+	    //
+	    //   Start the First Level Trigger simulation
+	    //
+	  
+	    if ( Lev0MT[ict] > 0 || Write_All_Images) {
+	      Lev1MT[ict]= Trigger_CT[ict]->FirstLevel();
 	    }
-	  
-	  Trigger.SetThreshold(fpixelthres);
-				 
-	  Trigger.Diskriminate() ; 
-
-	  //
-	  //   look if in all the signals in the trigger signal branch
-	  //   is a possible Trigger. Therefore we habe to diskriminate all
-	  //   the simulated analog signals (Method Diskriminate in class
-	  //   MTrigger). We look simultanously for the moments at which
-	  //   there are more than TRIGGER_MULTI pixels above the 
-	  //   CHANNEL_THRESHOLD. 
-	  //
-	  
-	  Lev0 = (Short_t) Trigger.ZeroLevel() ; 
-	  
-	  Lev1 = 0 ; 
-	  
-	  //
-	  //   Start the First Level Trigger simulation
-	  //
-	  
-	  if ( Lev0 > 0 || Write_All_Images) {
-	    Lev1= Trigger.FirstLevel();
+	    if (Lev1MT[ict]>0){
+	      ++ntrigger[ict];
+	    }
 	  }
-	  if (Lev1>0){
-	    ++ntrigger;
+
+	  Int_t NumImages = 0;
+	  Int_t CT_triggered=0;
+	  for(int ict=0;ict<ct_Number;ict++){
+	    if(NumImages==0 && Lev1MT[ict]>0)
+	      CT_triggered=ict;
+	    NumImages = (NumImages>=Lev1MT[ict])?NumImages:1;
+	    Lev0MT[ict]=1;
+	    if (Lev1MT[ict]==0 && Write_All_Images){ 
+	      NumImages=1;
+	      Lev0MT[ict]=0;
+	    }
 	  }
-	  Int_t NumImages = Lev1;
-	  Lev0=1;
-	  if (Lev1==0 && Write_All_Images){ 
-	    NumImages=1;
-	    Lev0=0;
-	  }
-
-	  for(Int_t ii=0;ii<NumImages;ii++){
-	    //  Loop over different level one triggers
-
+
+	  for(int ict=0;ict<ct_Number;ict++){
+	    for(Int_t ii=0;ii<NumImages;ii++){
+
+	      btrigger=1;
+	      
+	      //  Loop over different level one triggers
+	      
+	      //
+	      //   fill inside class fadc the member output
+	      //
+	      if(Lev1MT[ict]>0)
+		Fadc_CT[ict]->
+		  TriggeredFadc(Trigger_CT[ict]->GetFirstLevelTime(ii));
+	      else
+		Fadc_CT[ict]->
+		  TriggeredFadc(Trigger_CT[CT_triggered]->GetFirstLevelTime(ii));
+
+	      if(!flagstoring)
+		nstoredevents++;
+	      flagstoring = 1;
+
+	      if (Write_McTrig){
+		McTrig[ict]->SetFirstLevel (Lev1MT[ict]);
+		McTrig[ict]
+		  ->SetTime(Trigger_CT[ict]->GetFirstLevelTime(ii),ii+1);
+		Trigger_CT[ict]->GetMapDiskriminator(trigger_map);
+		McTrig[ict]->SetMapPixels(trigger_map,ii);
+	      }
+	      
+	      //  Fill Evt information
+	      
+	      if (Write_RawEvt){
+		
+		//
+		//  Fill the header of this event 
+		//
+		
+		EvtHeader[ict]
+		  ->FillHeader ( (UInt_t) (ntshow + nshow) ,  20 ) ; 
+	      
+		//   fill pixel information
+	      
+		if (Lev1MT[ict]){
+		  for(UInt_t i=0;
+		      i<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
+		      i++){
+		    for (j=0;j<FADC_SLICES;j++){
+		      fadcValues->AddAt(Fadc_CT[ict]->GetFadcSignal(i,j),j);
+		      fadcValuesLow->AddAt(Fadc_CT[ict]->GetFadcLowGainSignal(i,j),j);
+		    }
+		    EvtData[ict]->AddPixel(i,fadcValues,0);
+		    EvtData[ict]->AddPixel(i,fadcValuesLow,kTRUE);
+		  }
+		}
+	      }	
+	    }
 	    //
-	    //   fill inside class fadc the member output
-	    //
-	    fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
-
-	    nstoredevents++;
-
-	    if (Write_McTrig){
-	      McTrig[0]->SetFirstLevel ((ii+1)*Lev0);
-	      McTrig[0]->SetTime(Trigger.GetFirstLevelTime(ii),ii+1);
-	      Trigger.GetMapDiskriminator(trigger_map);
-	      McTrig[0]->SetMapPixels(trigger_map,ii);
+	    //    if a first level trigger occurred, then 
+	    //       1. do some other stuff (not implemented) 
+	    //       2. start the gui tool
+	    
+	    if(FADC_Scan){
+	      if ( Lev0MT[ict] > 0 ) {
+		  Fadc_CT[ict]->ShowSignal( McEvt, (Float_t) 60. ) ; 
+	      }
 	    }
-
-	    //  Fill Evt information
-
-	    if (Write_RawEvt){
-
-	      //
-	      //  Fill the header of this event 
-	      //
+	    
+	    if(Trigger_Scan){
+	      if ( Lev0MT[ict] > 0 ) {
+		Trigger_CT[ict]->ShowSignal(McEvt) ;
+	      }
+	    }
 	      
-	      EvtHeader[0]->FillHeader ( (UInt_t) (ntshow + nshow) ,  20 ) ; 
-	      
-	      //   fill pixel information
-	      
-	      if (Lev1){
-		for(UInt_t i=0;i<camgeom.GetNumPixels();i++){
-		  for (j=0;j<FADC_SLICES;j++){
-		    fadcValues->AddAt(fadc.GetFadcSignal(i,j),j);
-		    fadcValuesLow->AddAt(fadc.GetFadcLowGainSignal(i,j),j);
-		  }
-		  EvtData[0]->AddPixel(i,fadcValues,0);
-		  EvtData[0]->AddPixel(i,fadcValuesLow,kTRUE);
-		}
-	      }
-	    }	    
+	  } // end CT loop
+
+	  // If there is trigger in some telecope or we store all showers
+	  if(btrigger){  
 	    //
 	    //   fill the MMcEvt with all information  
@@ -1964,6 +2259,6 @@
 			     (UInt_t)(mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot()), 
 			     (UInt_t) ncph, 
-			     (UInt_t) inumphe,
-			     (UInt_t) inumphensb+inumphe,
+			     (UInt_t) inumphe_CT[0],
+			     (UInt_t) inumphensb+inumphe_CT[0],
 			     -1.0,
 			     -1.0,
@@ -2001,6 +2296,6 @@
 			     (UInt_t) (mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot()), 
 			     (UInt_t) ncph, 
-			     (UInt_t) inumphe,
-			     (UInt_t) inumphensb+inumphe,
+			     (UInt_t) inumphe_CT[0],
+			     (UInt_t) inumphensb+inumphe_CT[0],
 			     mcevth_2.get_ElecFraction(),
 			     mcevth_2.get_MuonFraction(),
@@ -2009,5 +2304,5 @@
 	    }
 	    //   We don not count photons out of the camera.	
-	    
+	      
 	    
 	    //
@@ -2016,28 +2311,13 @@
 	    
 	    EvtTree.Fill() ; 
-	    	  	  
-	    //
-	    //    if a first level trigger occurred, then 
-	    //       1. do some other stuff (not implemented) 
-	    //       2. start the gui tool
+	  }
 	    
-	    if(FADC_Scan){
-	      if ( Lev0 > 0 ) {
-		fadc.ShowSignal( McEvt, (Float_t) 60. ) ; 
-	      }
-	    }
-	    
-	    if(Trigger_Scan){
-	      if ( Lev0 > 0 ) {
-		Trigger.ShowSignal(McEvt) ;
-	      }
-	    }
-	    
-	    //    clear all
-	    if (Write_RawEvt) EvtHeader[0]->Clear() ;
-	    if (Write_RawEvt) EvtData[0]->DeletePixels();
-	    if (Write_McEvt) McEvt->Clear() ; 
+	  //    clear all
+	  for(int ict=0;ict<ct_Number;ict++){
+	    if (Write_RawEvt) EvtHeader[ict]->Clear() ;
+	    if (Write_RawEvt) EvtData[ict]->DeletePixels();
+	    if (Write_McTrig) McTrig[ict]->Clear() ;
 	  } 
-	  if (Write_McTrig) McTrig[0]->Clear() ;
+	  if (Write_McEvt) McEvt->Clear() ; 
 	}
 		
@@ -2064,5 +2344,6 @@
 	}
 	
-	for (int i=0; i<camgeom.GetNumPixels(); ++i) {
+	for (int i=0;
+	     i<((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels(); ++i) {
 	  printf("%d (%d): ", i, npixneig[i]);
 	  for (j=0; j<npixneig[i]; ++i) 
@@ -2073,8 +2354,12 @@
 #endif // __DEBUG__
 	  	  	  
+
+	// We search the maximum impact parameter fo the simualted showers
+	maxpimpact=maxpimpact<impactD?impactD:maxpimpact;
 	
 	// look for the next event
 	
-	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
+	for(int ict=0;ict<ct_Number;ict++)
+	  fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
 
       } // end while there is a next event
@@ -2088,5 +2373,7 @@
 	log(SIGNATURE, "End of this run with %d events . . .\n", nshow);
 	
-	fread( flag, SIZE_OF_FLAGS, 1, inputfile );
+	//fread( flag, SIZE_OF_FLAGS, 1, inputfile );
+	for(int ict=0;ict<ct_Number;ict++)
+	  fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
 	
 	if( isA( flag, FLAG_END_OF_FILE ) ){ // end of file
@@ -2094,13 +2381,13 @@
 	  still_in_loop  = FALSE;
 	  log(SIGNATURE, "Reading ASCII files at the end of the reflector file.  . .\n");
-	  read_ascii(inputfile, McConfigRunHeader);
-	  //break;
-	  if ((! Data_From_STDIN) && ( !feof(inputfile) )){
+	  for(int ict=0;ict<ct_Number;ict++){
+	    read_ascii(inputfile[ict], McConfigRunHeader[ict]);
+	  }
+	  if ((! Data_From_STDIN) && ( !feof(inputfile[0]) )){
 	    
 	    // we have concatenated input files.
 	    // get signature of the next part and check it.
 
-	    if((reflector_file_version=check_reflector_file( inputfile ))==FALSE){
-	      // exit(1);
+	    if((reflector_file_version=check_reflector_file( inputfile[0] ))==FALSE){
 	      log(SIGNATURE, "Next file is not recognised as a reflector file.\n");
 	      log(SIGNATURE, "Stopping ...\n");
@@ -2109,6 +2396,7 @@
 	    
 	  }
-	  
-	  fread( flag, SIZE_OF_FLAGS, 1, inputfile );
+
+	for(int ict=0;ict<ct_Number;ict++)	  
+	  fread( flag, SIZE_OF_FLAGS, 1, inputfile[ict] );
 
 	} // end if found end of file
@@ -2124,4 +2412,5 @@
   Float_t ftime;
   Float_t rnum;
+  Float_t  viewcone[2]={0,0};
 
   get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
@@ -2147,4 +2436,5 @@
       heights[i]=mcevth_2.get_HeightLev (i);
     rnum=mcevth_2.get_RunNumber();
+    mcevth_2.get_viewcone(&viewcone[0],&viewcone[1]);
   }
 
@@ -2166,5 +2456,5 @@
 		    Write_RawEvt,
 		    addElecNoise,
-		    CAMERA_PIXELS,
+		    ct_NPixels,
 		    (UInt_t)ntshow,
 		    (UInt_t)nstoredevents,
@@ -2185,4 +2475,5 @@
 		    shphimax,
 		    shphimin,
+		    maxpimpact,
 		    mcevth.get_CWaveLower(),
 		    mcevth.get_CWaveUpper(),
@@ -2206,5 +2497,5 @@
 		    Write_RawEvt,
 		    addElecNoise,
-		    CAMERA_PIXELS,
+		    ct_NPixels,
 		    (UInt_t)ntshow,
 		    (UInt_t)nstoredevents,
@@ -2225,4 +2516,5 @@
 		    shphimax,
 		    shphimin,
+		    maxpimpact,
 		    mcevth_2.get_CWaveLower(),
 		    mcevth_2.get_CWaveUpper(),
@@ -2257,5 +2549,5 @@
   Float_t constantNFL[4];
   mcrunh.get_constantNFL(&constantNFL[0]);
-
+  
   if(reflector_file_version>5)
     McCorsikaRunHeader->Fill(rnum,
@@ -2281,27 +2573,32 @@
 			     constantCATM,
 			     constantNFL,
+			     viewcone,
 			     mcrunh.get_wobble(),
 			     mcrunh.get_atmophere()
-			  );
-
-  //  Store qe for each PMT in output file
-  TArrayF qe_pmt;
-  TArrayF wav_pmt;
-  McConfigRunHeader->InitSizePMTs(ct_NPixels);
-  for(int i=0; i<ct_NPixels;i++){
-    McConfigRunHeader->AddPMT(i);
-    MGeomPMT &pmt  = McConfigRunHeader->GetPMT(i);
-    qe_pmt.Set(pointsQE,QE[i][1]);
-    wav_pmt.Set(pointsQE,QElambda);
-    pmt.SetArraySize(pointsQE);
-    pmt.SetPMTContent(i,wav_pmt,qe_pmt);
-  }
-
-  // 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->SetLightGuides(theta_lg,factor_lg);
+			     );
+  
+    //  Store qe for each PMT in output file
+    TArrayF qe_pmt;
+    TArrayF wav_pmt;
+
+  for(int ict=0;ict<ct_Number;ict++){
+    McConfigRunHeader[ict]->InitSizePMTs(ct_NPixels);
+    for(int i=0; i<(Int_t)((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();i++){
+      McConfigRunHeader[ict]->AddPMT(i);
+      MGeomPMT &pmt  = McConfigRunHeader[ict]->GetPMT(i);
+      qe_pmt.Set(pointsQE[ict],QE[ict][i][1]);
+      wav_pmt.Set(pointsQE[ict],QElambda);
+      pmt.SetArraySize(pointsQE[ict]);
+      pmt.SetPMTContent(i,wav_pmt,qe_pmt);
+    }
+    
+    // Store Light 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);
+  
+  }
 
   //  Fill the Header Tree with the current leaves of each branch
@@ -2336,7 +2633,9 @@
   }
   else{
-    log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n", 
-	 ((float)ntrigger) / ((float)ntshow) * 100.0, ntrigger, ntshow);
-    datafile<<"Fraction of triggers: "<<((float)ntrigger) / ((float)ntshow) * 100.0<<" ("<<ntrigger<<" out of "<<ntshow<<" )"<<endl;
+    for(int ict=0;ict<ct_Number;ict++){
+      log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n", 
+	   ((float)ntrigger[ict]) / ((float)ntshow) * 100.0, ntrigger[ict], ntshow);
+      datafile<<"Fraction of triggers: "<<((float)ntrigger[ict]) / ((float)ntshow) * 100.0<<" ("<<ntrigger[ict]<<" out of "<<ntshow<<" )"<<endl;
+    }
   }
 
@@ -2346,5 +2645,6 @@
   
   if( ! Data_From_STDIN ){
-    fclose( inputfile );
+    for(int ict=0;ict<ct_Number;ict++)
+      fclose( inputfile[ict] );
   }
   datafile.close();
@@ -2740,5 +3040,5 @@
 //!@{
 void 
-read_QE(char fname[256]){
+read_QE(char fname[256], int ict){
   ifstream qefile;
   char line[LINE_MAX_LENGTH];
@@ -2784,19 +3084,19 @@
       // get the number of datapoints 
 
-      sscanf(line, "%d", &pointsQE);
+      sscanf(line, "%d", &pointsQE[ict]);
       
       // allocate memory for the table of QEs
       
-      QE = new float ** [ct_NPixels];
+      QE[ict] = new float ** [ct_NPixels];
 
       for ( i=0; i<ct_NPixels; ++i ) {
-        QE[i] = new float * [2];
-        QE[i][0] = new float[pointsQE];
-        QE[i][1] = new float[pointsQE];
+        QE[ict][i] = new float * [2];
+        QE[ict][i][0] = new float[pointsQE[ict]];
+        QE[ict][i][1] = new float[pointsQE[ict]];
       }
       
-      QElambda = new float [pointsQE];
-
-      for ( i=0; i<pointsQE; ++i ) {
+      QElambda = new float [pointsQE[ict]];
+
+      for ( i=0; i<pointsQE[ict]; ++i ) {
         qefile.getline(line, LINE_MAX_LENGTH);
         sscanf(line, "%f", &QElambda[i]);
@@ -2814,7 +3114,7 @@
 
     if ( ((i-1) < ct_NPixels) && ((i-1) > -1) &&
-         ((j-1) < pointsQE)   && ((j-1) > -1) ) {
-      QE[i-1][0][j-1] = QElambda[j-1];
-      QE[i-1][1][j-1] = qe;
+         ((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;
     }
 
@@ -2825,7 +3125,7 @@
   }
 
-  if(icount/pointsQE < ct_NPixels){
+  if(icount/pointsQE[ict] < ct_NPixels){
     error( "read_QE", "The quantum efficiency file is faulty\n  (found only %d pixels instead of %d).\n",
-	   icount/pointsQE, ct_NPixels );
+	   icount/pointsQE[ict], ct_NPixels );
   }
 
@@ -2837,9 +3137,9 @@
 
   for(icount=0; icount< ct_NPixels; icount++){
-    for(i=0; i<pointsQE; i++){
-      if( QE[icount][0][i] < 100. || QE[icount][0][i] > 1000. ||
-	  QE[icount][1][i] < 0. || QE[icount][1][i] > 100.){
+    for(i=0; i<pointsQE[ict]; i++){
+      if( QE[ict][icount][0][i] < 100. || QE[ict][icount][0][i] > 1000. ||
+	  QE[ict][icount][1][i] < 0. || QE[ict][icount][1][i] > 100.){
 	error( "read_QE", "The quantum efficiency file is faulty\n  pixel %d, point %d  is % f, %f\n",
-	       icount, i, QE[icount][0][i], QE[icount][1][i] );
+	       icount, i, QE[ict][icount][0][i], QE[ict][icount][1][i] );
       }
     }
@@ -2928,5 +3228,5 @@
 read_pixels(struct camera *pcam)
 {
-  ifstream qefile;
+  /*  ifstream qefile;
   char line[LINE_MAX_LENGTH];
   int n, i, j, icount;
@@ -3103,5 +3403,5 @@
 
   log("read_pixels", "Done.\n");
-
+  */
 }
 //!@}
@@ -3148,7 +3448,17 @@
       fgets(line,500,sp);
       sscanf(line,"%i",&numref);
+      TArrayF wavarray(numref);
+      TArrayF refarray(numref);
+      fgets(line,500,sp);
       for(int i=0; i<numref;i++){
 	fgets(line,500,sp);
 	sscanf(line,"%f %f",&wav,&ref);
+	wavarray[i]=wav;
+	refarray[i]=ref;
+      }
+      for (int j=0; j<nummir;j++){
+	MGeomMirror &mirror = config->GetMirror(j);
+	mirror.SetArraySize(numref);
+	mirror.SetReflectivity(wavarray, refarray);
       }
       continue;
@@ -3533,5 +3843,5 @@
 #define sqrt3  1.732050807 // = sqrt(3.)
 
-int bpoint_is_in_pix(double dx, double dy, MGeomCamMagic *pgeomcam)
+int bpoint_is_in_pix(double dx, double dy, MGeomCam *pgeomcam)
 {
     /* return TRUE if point (dx, dy) is in pixel number ipixnum, else return FALSE (use camera coordinate system) */
@@ -3548,10 +3858,10 @@
     {
       MGeomPix &pixel = (*pgeomcam)[i];
-      const double b = pixel.GetR()/2;
-      const double a = pixel.GetR()/sqrt3;
+      const double b = pixel.GetD()/2;
+      const double a = pixel.GetD()/sqrt3;
       
       const double xx = dx - pixel.GetX();
       const double yy = dy - pixel.GetY();
-      
+
       if(((-b <= xx) && (xx <= 0.) && ((-sqrt13 * xx - a) <= yy) && (yy <= ( sqrt13 * xx + a))) ||
 	 ((0. <  xx) && (xx <= b ) && (( sqrt13 * xx - a) <= yy) && (yy <= (-sqrt13 * xx + a)))   ){
@@ -3761,5 +4071,5 @@
 
 int produce_phes( FILE *sp, // the input file
-		  class MGeomCamMagic *camgeom, // the camera layout
+		  class MGeomCam *camgeom, // the camera layout
 		  float minwl_nm, // the minimum accepted wavelength
 		  float maxwl_nm, // the maximum accepted wavelength
@@ -3770,5 +4080,6 @@
 		  int *incph,    // total number of cph read
 		  float *tmin_ns,   // minimum arrival time of all phes
-		  float *tmax_ns    // maximum arrival time of all phes
+		  float *tmax_ns,    // maximum arrival time of all phes
+		  int ict      // Telescope that is being analised to get the right QE.
                    ){
 
@@ -3866,9 +4177,9 @@
     // set pointer to the QE table of the relevant pixel
     	  
-    qept = (float **)QE[ipixnum];
+    qept = (float **)QE[ict][ipixnum];
 	  
     // check if wl is inside table; outside the table, QE is assumed to be zero
     
-    if((wl < qept[0][0]) || (wl > qept[0][pointsQE-1])){
+    if((wl < qept[0][0]) || (wl > qept[0][pointsQE[ict]-1])){
       continue;
 	    
@@ -3878,5 +4189,5 @@
 
     k = 1; // start at 1 because the condition was already tested for 0
-    while (k < pointsQE-1 && qept[0][k] < wl){
+    while (k < pointsQE[ict]-1 && qept[0][k] < wl){
       k++;
     }
@@ -3923,5 +4234,5 @@
     fadc->Fill(ipixnum,(t-*tmin_ns) , 
 	       trigger->FillShow(ipixnum,t-*tmin_ns),
-	       !((*camgeom)[ipixnum].GetR()>(*camgeom)[0].GetR()));
+	       !((*camgeom)[ipixnum].GetD()>(*camgeom)[0].GetD()));
     
     *itotnphe += 1;
@@ -3944,7 +4255,8 @@
 
 int produce_nsbrates( char *iname, // the starfield input file name
-		      MGeomCamMagic *camgeom, // camera layout
-		      float rate_phepns[iMAXNUMPIX][iNUMWAVEBANDS] // the product of this function:
+		      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
 		      ){
 
@@ -3952,7 +4264,7 @@
   int k, ii; // counters
 
-  MTrigger trigger(Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
+  MTrigger trigger(397,camgeom,Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
   MFadc flashadc;
-		      
+
   static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
 					    WAVEBANDBOUND2,
@@ -4100,5 +4412,6 @@
 			      &incph,
 			      &tmin_ns,
-			      &tmax_ns
+			      &tmax_ns,
+			      0
 			      ); // photons from starfield
 
@@ -4269,4 +4582,14 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.57  2003/07/17 18:02:46  blanch
+// Several new features introduced as well as fixed bugs
+//
+// 	- 1/100 events printed out
+// 	- Low gain implemented
+// 	- Different response for outer and inner pixels
+// 	- Some warnings removed
+// 	- pedestal and qe file from inpuit card
+// 	- Faster electronic simulation
+//
 // Revision 1.55  2003/02/12 12:22:10  blanch
 // *** empty log message ***
@@ -4436,179 +4759,4 @@
 // is used.
 //
-// Revision 1.12  2000/09/22 17:40:18  harald
-// Added a lot of changes done by oscar.
-//
-// Revision 1.11  2000/09/21 11:47:33  harald
-// Oscar found some smaller errors in the calculation of the pixel shape and
-// corrected it.
-//
-// $Log: not supported by cvs2svn $
-// Revision 1.55  2003/02/12 12:22:10  blanch
-// *** empty log message ***
-//
-// Revision 1.54  2003/02/12 11:55:01  blanch
-// *** empty log message ***
-//
-// Revision 1.53  2003/01/23 18:35:21  blanch
-// *** empty log message ***
-//
-// Revision 1.52  2003/01/20 17:19:20  blanch
-// It fills the MMcCorsikaRun.
-//
-// Revision 1.51  2003/01/14 13:40:17  blanch
-// MRawRunHeader::fNumEvents has been filled with the number of events in
-// this file.
-// Problems in fImpact computation have been solved.
-// Option to set a dc value to rise the discriminator threshold has been added.
-//
-// Revision 1.50  2003/01/07 16:33:31  blanch
-// Star Field Rotation has been implemented by Raul Orduna. Now there is a
-// rotation for each shower. It is done by a non enter pixel rotation assuming
-// a circular symetry of the camera. It is not exact but it is accurate enough and
-// much faster.
-//
-// Revision 1.49  2002/12/13 10:04:07  blanch
-// *** empty log message ***
-//
-// Revision 1.48  2002/12/12 17:40:50  blanch
-// *** empty log message ***
-//
-// Revision 1.47  2002/12/10 17:19:31  blanch
-// *** empty log message ***
-//
-// Revision 1.46  2002/11/08 17:51:00  blanch
-// *** empty log message ***
-//
-// Revision 1.45  2002/10/29 17:15:27  blanch
-// Reading several reflector versions.
-//
-// Revision 1.44  2002/10/18 16:53:03  blanch
-// Modification to read several reflector files.
-//
-// Revision 1.43  2002/09/13 10:53:39  blanch
-// Minor change to remove some undisired comments.
-//
-// Revision 1.42  2002/09/09 16:00:49  blanch
-// Statement has been included to avoid writting to disk MParContainer and MArray.
-// It has also been added the effect of the WC, the actual values must be added,
-// once they are measured.
-//
-// Revision 1.41  2002/09/04 09:57:42  blanch
-// Modifications done to use MGeomCam from MARS.
-//
-// Revision 1.40  2002/07/16 16:15:22  blanch
-// A first implementation of the Star field rotation has been done.
-//
-// Revision 1.39  2002/04/27 10:48:39  blanch
-// Some unused varibles have been removed.
-//
-// Revision 1.38  2002/03/18 18:44:29  blanch
-// Small modificatin to set the electronic Noise in the MMcTrigHeader class.
-//
-// Revision 1.37  2002/03/18 16:42:20  blanch
-// The data member fProductionSite of the MMcRunHeader has been set to 0,
-// which means that the prodution site is unknown.
-//
-// Revision 1.35  2002/03/15 15:15:52  blanch
-// Several modifications to simulate the actual trigger zone.
-//
-// Revision 1.34  2002/03/13 18:13:56  blanch
-// Some changes to fill correctly the new format of MMcRunHeader.
-//
-// Revision 1.33  2002/03/04 17:21:48  blanch
-// Small and not important changes.
-//
-// Revision 1.32  2002/02/28 15:04:52  blanch
-// A small back has been solved. Before, while not using the option
-// writte_all_images, not all triggered showers were stored. Now it is solved.
-// For that it is importatn taht the less restrictive trigger option is
-// checked first.
-// A new facility has been introduced and now one can choose the step size in
-// trigger loop mode for the discriminator threshold.
-// The close-compact topology for less than 3 pixels does not make sense. Before
-// the program was ignoring that, now it switch to simple neighbour condition.
-//
-// Revision 1.31  2002/01/18 17:41:02  blanch
-// The option of adding noise to all pixels or to not adding the noise
-// has been added.
-// We removed the pixels larger than 577. When there were more than one
-// trigger in one shower, the pixel number was increasing. Now it is
-// flagged by the variable MMcTrig::fFirstLvlTrig.
-//
-// Revision 1.30  2001/11/27 09:49:54  blanch
-// Fixing bug which was treating wrongly the extension of star photons.
-//
-// Revision 1.29  2001/11/14 17:38:23  blanch
-// Sveral changes have been done:
-// 	- bpoint_in_in_pixel has been dodified to speed up the program
-// 	- Some minor changes have been done to remove warnings
-// 	- buffer size and split version of the Branches have been removed
-// 	- Some modifications were needed to adat the program to the new
-// 	  MRawEvtData::DeletePixels
-//
-// Revision 1.28  2001/10/26 16:31:45  blanch
-// Removing several warnings.
-//
-// Revision 1.27  2001/09/05 10:04:33  blanch
-// *** empty log message ***
-//
-// Revision 1.26  2001/07/19 09:29:53  blanch
-// Different threshold for each pixel can be used.
-//
-// Revision 1.25  2001/05/08 08:07:54  blanch
-// New numbering for branches from different trigger conditions has been
-// implemented. Now, they are calles: ClassName;1., ClasNema;2., ...
-// The MontaCarlo Headers (MMcTrigHeader and MMcFadcHeader) have been move to
-// the RunHeaders tree. Also the MRawRunHeader is thera with some of its
-// information already filled.
-//
-// Revision 1.24  2001/04/06 16:48:09  magicsol
-// New camera version able to read the new format of the reflector output:
-// reflector 0.4
-//
-// Revision 1.23  2001/03/28 16:13:41  blanch
-// While feeling the MMcFadeHeader branch the boolean conditoin was wrong. It has
-// been solved.
-//
-// Revision 1.22  2001/03/20 18:52:43  blanch
-// *** empty log message ***
-//
-// Revision 1.21  2001/03/19 19:53:03  blanch
-// Some print outs have been removed.
-//
-// Revision 1.20  2001/03/19 19:30:06  magicsol
-// Minor changes have been done to improve the FADC pedestals treatment.
-//
-// Revision 1.19  2001/03/05 11:14:41  magicsol
-// I changed the position of readinf a parameter. It is a minnor change.
-//
-// Revision 1.18  2001/03/05 10:36:52  blanch
-// A branch with information about the FADC simulation (MMcFadcHeader) is writen
-// in the McHeaders tree of the aoutput root file.
-// The NSB contribution is only applied if the the number of phe form the shower
-// are above a given number.
-//
-// Revision 1.17  2001/02/23 11:05:57  magicsol
-// Small changes due to slightly different output format and the introduction of
-// pedesals for teh FADC.
-//
-// Revision 1.16  2001/01/18 18:44:40  magicsol
-// *** empty log message ***
-//
-// Revision 1.15  2001/01/17 09:32:27  harald
-// The changes are neccessary to have the same name for trees in MC and in
-// data. So now there should be no differences in MC data and real data from
-// FADC system.
-//
-// Revision 1.14  2001/01/15 12:33:34  magicsol
-// Some modifications have been done to use the new (Dec'2000) Raw data format.
-// There are also some minnors modifications to adapt the improvements in the
-// MTrigger class (overlaping time and trigger cells).
-//
-// Revision 1.13  2000/10/25 08:14:23  magicsol
-// The routine that produce poisson random numbers to decide how many phe
-// form NSB are emmited in each pixel has been replaced. Now a ROOT routine
-// is used.
-//
 // Revision 1.10  2000/07/04 14:10:20  MagicSol
 // Some changes have been done in the root output file. The RawEvt tree is only