Changeset 2334 for trunk/MagicSoft/Simulation/Detector

Timestamp:

09/15/03 10:59:53 (21 years ago)

Author:

blanch

Message:

The concept of the camera prgoram has not changed but this version has
quite a lot of changes to allow several Camera geometries as well as
multitelescope.

It is suposed to be the first working code for camera 0.7.

File:

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

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

@@ -21 +21 @@
 //
 // $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 +58 @@
 
 #include "TArrayC.h"
-
+#include "TObjArray.h"
+                                                         
 #include "MTrigger.hxx"
 #include "MFadc.hxx"
@@ -75 +76 @@
 #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 +196 @@
 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 +258 @@
 
 //@: 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 +314 @@
 
 //@: 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 +411 @@
   //@'
 
-  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 +427 @@
   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 +442 @@
 
   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 +459 @@
   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 +472 @@
   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 +497 @@
   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 +532 @@
   Float_t soffphi = 0.0;
   UInt_t corsika = 5200 ;
+  Float_t maxpimpact = 0.0;
 
   // Star Field Rotation variables
@@ -522 +539 @@
   Float_t rho  , C3 , C2 , C1; 
 
-  MGeomCamMagic camgeom; // structure holding the camera definition
-
-  ct_NPixels=camgeom.GetNumPixels();
-
   //!@' @#### Definition of variables for |getopt()|.
   //@'
@@ -534 +547 @@
   //@'
 
+  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 +616 @@
   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 +705 @@
                        &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 +715 @@
     (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 +770 @@
       "%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 +794 @@
         "%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 +803 @@
         "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 +890 @@
 
   Int_t Lev0, Lev1; 
+  Int_t Lev0MT[100], Lev1MT[100]; 
   
   fadcValues    =  new TArrayC(FADC_SLICES);
@@ -804 +897 @@
   // number of pixels for parameters
     
-  anaPixels = get_ana_pixels();
-  anaPixels = (anaPixels == -1) ? camgeom.GetNumPixels() : anaPixels;
-
   // Switched off writing TObject
 
@@ -818 +908 @@
   // 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 +942 @@
   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 +966 @@
   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 +995 @@
   /////////////////////////////////////////////////////////////////////
 
-  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 +1029 @@
 
   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 +1041 @@
 
   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 +1071 @@
                     &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 +1106 @@
   }  
 
-  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 +1132 @@
   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 +1147 @@
     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 +1167 @@
           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 +1183 @@
     }
   }
-
+  
   //  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 +1219 @@
           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 +1227 @@
   }
 
-  //  Fill the Header Tree with the current leaves of each branch
-  // HeaderTree.Fill() ;
-            
 
   //      create a Tree for the Event data stream 
@@ -1085 +1237 @@
   }
 
-  if(!Trigger_Loop){
+  if(!Trigger_Loop && ct_Number==1){
     
+    HeaderTree.Branch("MMcConfigRunHeader","MMcConfigRunHeader",
+                      &McConfigRunHeader[0]);
+
     if (Write_RawEvt){
       EvtTree.Branch("MRawEvtHeader","MRawEvtHeader", 
@@ -1099 +1254 @@
   }
   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 +1279 @@
   }  
 
-  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 +1289 @@
                      &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 +1300 @@
   }
 
+
   TApplication theAppTrigger("App", &argc, argv);
 
@@ -1171 +1339 @@
     //
     
+    // 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 +1350 @@
     }
   
-    // 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 +1413 @@
   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 +1452 @@
   // allocate space for PMTs numbers of pixels
   
-  fnpix = new float [ camgeom.GetNumPixels() ];
+  fnpix = new float [ct_NPixels];
 
      
@@ -1252 +1459 @@
 
   // 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 +1468 @@
   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 +1486 @@
     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 +1535 @@
         // 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 +1680 @@
         }
 
-        // 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 +1687 @@
           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 +1743 @@
 
             // 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 +1833 @@
         //   
 
-        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 +1855 @@
           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 +1886 @@
             
             //  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 +1906 @@
                   // 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 +1919 @@
                   //
                   if(Lev0!=0)
-                      Lev1=Trigger.FirstLevel();
+                      Lev1=Trigger_CT[0]->FirstLevel();
                   if(Lev1>0) {
                     btrigger= 1;
@@ -1694 +1936 @@
                       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 +1944 @@
                     //
 
-                    fadc.TriggeredFadc(Trigger.GetFirstLevelTime(ii));
+                    Fadc_CT[0]->TriggeredFadc(Trigger_CT[0]->GetFirstLevelTime(ii));
             
                     if( Write_RawEvt ){
@@ -1712 +1954 @@
                       //   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 +2020 @@
                              (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 +2057 @@
                              (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 +2069 @@
             //  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 +2086 @@
         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 +2259 @@
                              (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 +2296 @@
                              (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 +2304 @@
             }
             //   We don not count photons out of the camera.    
-            
+              
             
             //
@@ -2016 +2311 @@
             
             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 +2344 @@
         }
         
-        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 +2354 @@
 #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 +2373 @@
         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 +2381 @@
           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 +2396 @@
             
           }
-          
-          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 +2412 @@
   Float_t ftime;
   Float_t rnum;
+  Float_t  viewcone[2]={0,0};
 
   get_starfield_center(&sfRaH,&sfRaM,&sfRaS,&sfDeD,&sfDeM,&sfDeS);
@@ -2147 +2436 @@
       heights[i]=mcevth_2.get_HeightLev (i);
     rnum=mcevth_2.get_RunNumber();
+    mcevth_2.get_viewcone(&viewcone[0],&viewcone[1]);
   }
 
@@ -2166 +2456 @@
                     Write_RawEvt,
                     addElecNoise,
-                    CAMERA_PIXELS,
+                    ct_NPixels,
                     (UInt_t)ntshow,
                     (UInt_t)nstoredevents,
@@ -2185 +2475 @@
                     shphimax,
                     shphimin,
+                    maxpimpact,
                     mcevth.get_CWaveLower(),
                     mcevth.get_CWaveUpper(),
@@ -2206 +2497 @@
                     Write_RawEvt,
                     addElecNoise,
-                    CAMERA_PIXELS,
+                    ct_NPixels,
                     (UInt_t)ntshow,
                     (UInt_t)nstoredevents,
@@ -2225 +2516 @@
                     shphimax,
                     shphimin,
+                    maxpimpact,
                     mcevth_2.get_CWaveLower(),
                     mcevth_2.get_CWaveUpper(),
@@ -2257 +2549 @@
   Float_t constantNFL[4];
   mcrunh.get_constantNFL(&constantNFL[0]);
-
+  
   if(reflector_file_version>5)
     McCorsikaRunHeader->Fill(rnum,
@@ -2281 +2573 @@
                              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 +2633 @@
   }
   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 +2645 @@
   
   if( ! Data_From_STDIN ){
-    fclose( inputfile );
+    for(int ict=0;ict<ct_Number;ict++)
+      fclose( inputfile[ict] );
   }
   datafile.close();
@@ -2740 +3040 @@
 //!@{
 void 
-read_QE(char fname[256]){
+read_QE(char fname[256], int ict){
   ifstream qefile;
   char line[LINE_MAX_LENGTH];
@@ -2784 +3084 @@
       // 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 +3114 @@
 
     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 +3125 @@
   }
 
-  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 +3137 @@
 
   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 +3228 @@
 read_pixels(struct camera *pcam)
 {
-  ifstream qefile;
+  /*  ifstream qefile;
   char line[LINE_MAX_LENGTH];
   int n, i, j, icount;
@@ -3103 +3403 @@
 
   log("read_pixels", "Done.\n");
-
+  */
 }
 //!@}
@@ -3148 +3448 @@
       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 +3843 @@
 #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 +3858 @@
     {
       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 +4071 @@
 
 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 +4080 @@
                   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 +4177 @@
     // 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 +4189 @@
 
     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 +4234 @@
     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 +4255 @@
 
 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 +4264 @@
   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 +4412 @@
                               &incph,
                               &tmin_ns,
-                              &tmax_ns
+                              &tmax_ns,
+                              0
                               ); // photons from starfield
 
@@ -4269 +4582 @@
 //
 // $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 +4759 @@
 // 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