Index: trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx
===================================================================
--- trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 2418)
+++ trunk/MagicSoft/Simulation/Detector/Camera/camera.cxx	(revision 2427)
@@ -21,7 +21,7 @@
 //
 // $RCSfile: camera.cxx,v $
-// $Revision: 1.64 $
+// $Revision: 1.65 $
 // $Author: blanch $ 
-// $Date: 2003-10-21 07:42:50 $
+// $Date: 2003-10-26 19:43:00 $
 //
 ////////////////////////////////////////////////////////////////////////
@@ -165,7 +165,4 @@
 static int Write_All_Images = FALSE;
 
-//@: flag: TRUE: write all data to output; FALSE: only triggered showers
-static int Write_All_Data = FALSE;
-
 static int Write_McEvt  = TRUE;
 static int Write_McTrig = TRUE;
@@ -342,19 +339,19 @@
 
   MCRunHeader mcrunh;                //@< Run Header class (MC)
-  MCEventHeader mcevth;       //@< Event Header class (MC)
-  MCEventHeader_2 mcevth_2;   //@< Event Header class (MC) for reflector > 0.6
+  MCEventHeader mcevth[100];       //@< Event Header class (MC)
+  MCEventHeader_2 mcevth_2[100];   //@< Event Header class (MC) for reflector > 0.6
   MCCphoton cphoton;          //@< Cherenkov Photon class (MC)
 
   int inumphe;                //@< number of photoelectrons in an event from showers
   int inumphe_CT[100];         //@< number of photoelectrons in an event from showers
-  float inumphensb;             //@< number of photoelectrons in an event from nsb
+  float inumphensb[100];             //@< number of photoelectrons in an event from nsb
 
   float arrtmin_ns;           //@ arrival time of the first photoelectron
   float arrtmax_ns;           //@ arrival time of the last photoelectron
 
-  float thetaCT, phiCT;       //@< parameters of a given shower
+  float thetaCT[100], phiCT[100];       //@< parameters of a given shower
   float thetashw, phishw;     //@< parameters of a given shower
   float coreX, coreY;         //@< core position
-  float impactD;              //@< impact parameter
+  float impactD[100];              //@< impact parameter
   float l1, m1, n1;           //@< auxiliary variables
   float l2, m2, n2;           //@< auxiliary variables
@@ -364,9 +361,8 @@
   int nshow=0;                //@< partial number of shower in a given run
   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
-
-  int i, j, k;                //@< simple counters
+  int ncph[100];            //@< partial number of photons in a given run
+  int ntcph[100];                //@< total number of photons
+
+  int ibr, j, k;                //@< simple counters
 
   int addElecNoise;           //@< Will we add ElecNoise?
@@ -391,8 +387,9 @@
   Byte_t trigger_map[((Int_t)(CAMERA_PIXELS/8))+1]; //@< Pixels on when the
                                                     //@< camera triggers
-  Float_t fadc_elecnoise[CAMERA_PIXELS];  //@< Electronic niose for each pixel
+  Float_t fadc_elecnoise[CAMERA_PIXELS];  //@< Electronic noise for each pixel
+  Float_t fadc_diginoise[CAMERA_PIXELS];  //@< Digital noise for each pixel
   Float_t fadc_pedestals[CAMERA_PIXELS];  //@< array for fadc pedestals values
   Float_t fadc_sigma[CAMERA_PIXELS];      //@< array for fadc pedestals sigma
-  Float_t fadc_sigma_low[CAMERA_PIXELS];      //@< array for fadc pedestals sigma
+  Float_t fadc_sigma_low[CAMERA_PIXELS];  //@< array for fadc pedestals sigma
 
   float ext[iNUMWAVEBANDS] = { //@< average atmospheric extinction in each waveband
@@ -456,8 +453,8 @@
 
   qThreshold = new float *[100];
-  for(i=0;i<100;i++)
+  for(int i=0;i<100;i++)
     qThreshold[i] = new float [CAMERA_PIXELS]; 
 
-  for(i=0;i<iMAXNUMPIX;i++){
+  for(int i=0;i<iMAXNUMPIX;i++){
     for(int ict=0;ict<100;ict++){      
       nsb_phepns[ict][i]=0;
@@ -467,4 +464,6 @@
       nsbrate_phepns[i][j]=0.0;    //@< Starlight rates initialised at 0.0
   }
+  for(int i=0;i<100;i++)
+    ntcph[i]=0;
   /*!@'
 
@@ -529,11 +528,11 @@
   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;
+  for(int ict=0;ict<ct_Number;ict++){
+    ntrigger[ict]=0;
+    if(ct_Geometry>=ict*10){
+      GeometryCamera[ict]=int(ct_Geometry/pow(10.0,ict))%10;
     }
     else
-      GeometryCamera[i]=0;
+      GeometryCamera[ict]=0;
   }
 
@@ -603,8 +602,4 @@
   Write_All_Images = get_write_all_events();
 
-  // write all data (i.e., ph.e.s in pixels)
-
-  Write_All_Data = get_write_all_data();
-
   Write_McEvt  = get_write_McEvt()  ; 
   Write_McTrig = get_write_McTrig() ; 
@@ -752,11 +747,11 @@
   
   log(SIGNATURE,
-      "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %3.2f(%s) %3.2f(%s) %s\n",
+      "%s:\n\t%20s: %s\n\t%20s: %s\n\t%20s: %s\n\t%20s: %3.2f(%s) %3.2f+%3.2f(%s) %s\n",
       "Flags",
       "Data_From_STDIN",   ONoff(Data_From_STDIN),  
       "Write_All_Events",  ONoff(Write_All_Images), 
-      "Write_All_Data",    ONoff(Write_All_Data),
       "Rotate Starfield",  ONoff(Starfield_rotate),
-      "Electronic Noise",  Trigger_noise,"trigger", FADC_noise,"fadc",ONoff(addElecNoise)
+      "Electronic Noise",  Trigger_noise,"trigger", 
+      FADC_noise,DIGITAL_noise,"fadc",ONoff(addElecNoise)
       );
 
@@ -815,5 +810,5 @@
 
     log(SIGNATURE, "There are some showers to skip:\n");
-    for (i=0; i<nSkip; ++i)
+    for (int i=0; i<nSkip; ++i)
       log(SIGNATURE, "\tshower # %d\n", Skip[i]);
   }
@@ -890,5 +885,5 @@
     McTrig = new MMcTrig * [numberBranches];
   
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       McTrig[i] = new MMcTrig();
     }
@@ -896,5 +891,5 @@
     HeaderTrig = new MMcTrigHeader * [numberBranches];
   
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       HeaderTrig[i] = new MMcTrigHeader();
     }
@@ -905,5 +900,5 @@
     HeaderFadc = new MMcFadcHeader * [numberBranches];
   
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       HeaderFadc[i] = new MMcFadcHeader();
     }
@@ -934,5 +929,5 @@
   Float_t input_pedestals[ct_NPixels];
 
-  for(i=0;i<ct_NPixels;i++)
+  for(int i=0;i<ct_NPixels;i++)
     input_pedestals[i]=get_FADC_pedestal();
   for (int ict=0; ict<ct_Number;ict++){
@@ -961,5 +956,5 @@
   MMcConfigRunHeader  **McConfigRunHeader = NULL;
   McConfigRunHeader = new MMcConfigRunHeader * [numberBranches];
-  for (i=0;i<numberBranches;i++) {
+  for (int i=0;i<numberBranches;i++) {
     McConfigRunHeader[i] = new MMcConfigRunHeader();
   }
@@ -972,5 +967,5 @@
     EvtHeader = new MRawEvtHeader * [numberBranches]; 
 
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       EvtHeader[i] = new MRawEvtHeader();
     }
@@ -984,5 +979,5 @@
     EvtData = new MRawEvtData * [numberBranches]; 
 
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       EvtData[i] = new MRawEvtData();
       EvtData[i]->Init(RunHeader);     //  We need thr RunHeader to read
@@ -991,6 +986,17 @@
   }
 
-  MMcEvt  *McEvt = new MMcEvt (); 
-
+  MMcEvt  **McEvt = NULL;
+
+  if (Write_McEvt) {
+    if (ct_Number>1){
+      McEvt = new MMcEvt *[numberBranches]; 
+      for (int i=0;i<numberBranches;i++)
+	McEvt[i]=new MMcEvt();
+    }
+    else {
+      McEvt = new MMcEvt *[1]; 
+      McEvt[0]=new MMcEvt();
+    }
+  }
   //  
   // initalize a temporal ROOT file 
@@ -1048,12 +1054,13 @@
 
   if ((Trigger_Loop || ct_Number>1) && Write_McTrig){
-    for(char branchname[20],i=0;i<numberBranches;i++){
+    ibr=0;
+    for(char branchname[20];ibr<numberBranches;ibr++){
       
-      sprintf(help,"%i",i+1);
+      sprintf(help,"%i",ibr+1);
       strcpy (branchname, "MMcTrigHeader;");
       strcat (branchname, & help[0]);
       strcat (branchname, ".");
       HeaderTree.Branch(branchname,"MMcTrigHeader", 
-			&HeaderTrig[i]);
+			&HeaderTrig[ibr]);
     }
   }  
@@ -1065,12 +1072,13 @@
   }
   if ((Trigger_Loop || ct_Number>1) && Write_McFADC){
-    for(char branchname[20],i=0;i<numberBranches;i++){
+    ibr=0;
+    for(char branchname[20];ibr<numberBranches;ibr++){
       
-      sprintf(help,"%i",i+1);
+      sprintf(help,"%i",ibr+1);
       strcpy (branchname, "MMcFadcHeader;");
       strcat (branchname, & help[0]);
       strcat (branchname, ".");
       HeaderTree.Branch(branchname,"MMcFadcHeader", 
-			&HeaderFadc[i]);
+			&HeaderFadc[ibr]);
     }
   }  
@@ -1139,6 +1147,7 @@
   //  Fill branches for MMcFadcHeader
 
-  for(i=0;i<ct_NPixels;i++){
+  for(int i=0;i<ct_NPixels;i++){
       fadc_elecnoise[i]=FADC_noise;
+      fadc_diginoise[i]=DIGITAL_noise;
   }
 
@@ -1158,5 +1167,5 @@
     HeaderFadc[0]->SetPedestal(&fadc_pedestals[0],
 			       ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
-    HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0],
+    HeaderFadc[0]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
 				((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels());
   }
@@ -1174,5 +1183,6 @@
       HeaderFadc[i]->SetLow2High(FADC_high2low);
       HeaderFadc[i]->SetPedestal(&fadc_pedestals[0],((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
-      HeaderFadc[i]->SetElecNoise(&fadc_elecnoise[0],((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
+      HeaderFadc[i]->SetElecNoise(&fadc_elecnoise[0], &fadc_diginoise[0],
+				  ((MGeomCam*)(camgeom.UncheckedAt(i)))->GetNumPixels());
     }
   }  
@@ -1194,5 +1204,6 @@
 	  HeaderFadc[iconcount]->SetLow2High(FADC_high2low);
 	  HeaderFadc[iconcount]->SetPedestal(&fadc_pedestals[0], ct_NPixels);
-	  HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0], ct_NPixels);
+	  HeaderFadc[iconcount]->SetElecNoise(&fadc_elecnoise[0],
+					      &fadc_diginoise[0],ct_NPixels);
 	  iconcount++;
 	}
@@ -1212,8 +1223,21 @@
   TTree EvtTree("Events","Normal Triggered Events");
 
-  if (Write_McEvt){
+  if (Write_McEvt && ct_Number==1){
 
     EvtTree.Branch("MMcEvt","MMcEvt", 
-		   &McEvt);
+		   &McEvt[0]);
+  }
+
+  if (Write_McEvt && ct_Number!=1){
+    ibr=0;
+    for(char branchname[10];ibr<numberBranches;ibr++){
+      
+      sprintf(help,"%i",ibr+1);
+      strcpy (branchname, "MMcEvt;");
+      strcat (branchname, & help[0]);
+      strcat (branchname, ".");
+      EvtTree.Branch(branchname,"MMcEvt", 
+		     &McEvt[ibr]);
+    }
   }
 
@@ -1234,12 +1258,13 @@
 
     if (Write_McTrig){
-      for(char branchname[10],i=0;i<numberBranches;i++){
+      ibr=0;
+      for(char branchname[10];ibr<numberBranches;ibr++){
       
-	sprintf(help,"%i",i+1);
+	sprintf(help,"%i",ibr+1);
 	strcpy (branchname, "MMcTrig;");
 	strcat (branchname, & help[0]);
 	strcat (branchname, ".");
 	EvtTree.Branch(branchname,"MMcTrig", 
-		       &McTrig[i]);
+		       &McTrig[ibr]);
       }
     }
@@ -1247,21 +1272,23 @@
 
   if ((Trigger_Loop || ct_Number>1) && Write_RawEvt){
-    for(char branchname[15],i=0;i<numberBranches;i++){
+    ibr=0;
+    for(char branchname[15];ibr<numberBranches;ibr++){
       
-      sprintf(help,"%i",i+1);
+      sprintf(help,"%i",ibr+1);
       strcpy (branchname, "MRawEvtHeader;");
       strcat (branchname, & help[0]);
       strcat (branchname, ".");
       EvtTree.Branch(branchname,"MRawEvtHeader", 
-		     &EvtHeader[i]);
-    }
-    for(char branchname[15],i=0;i<numberBranches;i++){
+		     &EvtHeader[ibr]);
+    }
+    ibr=0;
+    for(char branchname[15];ibr<numberBranches;ibr++){
       
-      sprintf(help,"%i",i+1);
+      sprintf(help,"%i",ibr+1);
       strcpy (branchname, "MRawEvtData;");
       strcat (branchname, & help[0]);
       strcat (branchname, ".");
       EvtTree.Branch(branchname,"MRawEvtData", 
-		     &EvtData[i]);
+		     &EvtData[ibr]);
     }
   }
@@ -1307,4 +1334,6 @@
     
     // FIXME --- star NSB different for each camera?
+    log(SIGNATURE,"Produce NSB rates from Star Field");
+
     k = produce_nsbrates( starfieldname,
 			  ((MGeomCam*)(camgeom.UncheckedAt(0))),
@@ -1352,5 +1381,5 @@
 	const Float_t size=
 	  (*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD();
-	diffnsb_phepns[ict][ui] = (Int_t(meanNSB*factorNSB_ct*100*size*size+0.5))/(100.0*size*size) * size*size*factorqe_NSB[ict];
+	diffnsb_phepns[ict][ui] = (Int_t(meanNSB*factorNSB_ct*100*size*size+0.5))/(100.0)*factorqe_NSB[ict];
       }
     }
@@ -1375,6 +1404,6 @@
   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   //
-  //  Now a empty event with the conditin in which the camera is run
-  //  is simulated. IN this way one gets an estimation of the 
+  //  Now 500 empty event with the conditin in which the camera is run
+  //  are simulated. IN this way one gets an estimation of the 
   //  sigma for the pedestal of each FADC channel.
   //  This sigma computtaion is done assuming any noise taht effects 
@@ -1385,36 +1414,44 @@
 
   for(int ict=0;ict<ct_Number;ict++){
-    Fadc_CT[ict]->Reset() ; 
-    if (addElecNoise){	
-      Fadc_CT[ict]->ElecNoise(FADC_noise) ;
-    }
-    if(simulateNSB){
+    for(UInt_t ui=0;
+	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
+	ui++){
+      fadc_sigma[ui]=0.0;
+      fadc_sigma_low[ui]=0.0;
+    }
+    for(int ie=0;ie<500;ie++){
+      Fadc_CT[ict]->Reset() ; 
+      if (addElecNoise){	
+	Fadc_CT[ict]->ElecNoise(FADC_noise) ;
+      }
+      if(simulateNSB){
+	for(UInt_t ui=0;
+	    ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
+	    ui++){
+	  if(nsb_phepns[ict][ui]>0.0){
+	    if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD()>(*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD()){
+	      k=lons_outer.GetResponse(nsb_phepns[ict][ui],0.01,
+				       & nsb_trigresp[0],
+				       & nsb_fadcresp[0]);
+	    }
+	    else{
+	      k=lons.GetResponse(nsb_phepns[ict][ui],0.01,
+				 & nsb_trigresp[0],& nsb_fadcresp[0]);
+	    }
+	    if(k==0){
+	      cout << "Exiting.\n";
+	      exit(1);
+	    }
+	    Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
+	  }
+	}
+      }
+      Fadc_CT[ict]->Pedestals();
       for(UInt_t ui=0;
 	  ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
 	  ui++){
-	if(nsb_phepns[ict][ui]>0.0){
-	  if((*((MGeomCam*)(camgeom.UncheckedAt(ict))))[ui].GetD()>(*((MGeomCam*)(camgeom.UncheckedAt(ict))))[0].GetD()){
-	    k=lons_outer.GetResponse(nsb_phepns[ict][ui],0.01,
-				     & nsb_trigresp[0],
-				     & nsb_fadcresp[0]);
-	  }
-	  else{
-	    k=lons.GetResponse(nsb_phepns[ict][ui],0.01,
-			       & nsb_trigresp[0],& nsb_fadcresp[0]);
-	  }
-	  if(k==0){
-	    cout << "Exiting.\n";
-	    exit(1);
-	  }
-	  Fadc_CT[ict]->AddSignal(ui,nsb_fadcresp);
-	}
+	fadc_sigma[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,1))/500.0;
+	fadc_sigma_low[ui]+=(Fadc_CT[ict]->GetPedestalNoise(ui,0))/500.0;
       }
-    }
-    Fadc_CT[ict]->Pedestals();
-    for(UInt_t ui=0;
-	ui<((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();
-	ui++){
-      fadc_sigma[ui]=Fadc_CT[ict]->GetPedestalNoise(ui,1);
-      fadc_sigma_low[ui]=Fadc_CT[ict]->GetPedestalNoise(ui,0);
     }
     HeaderFadc[ict]->SetPedestalSigma(&fadc_sigma_low[0],&fadc_sigma[0],
@@ -1554,9 +1591,11 @@
 	if (reflector_file_version<6){
 	  for(int ict=0;ict<ct_Number;ict++)
-	    fread( (char*)&mcevth, mcevth.mysize(), 1, inputfile[ict] );
+	    fread( (char*)&mcevth[ict], mcevth[ict].mysize(),
+		   1, inputfile[ict] );
 	}
 	else{
 	  for(int ict=0;ict<ct_Number;ict++)
-	    fread( (char*)&mcevth_2, mcevth_2.mysize(), 1, inputfile[ict] );
+	    fread( (char*)&mcevth_2[ict], mcevth_2[ict].mysize(),
+		   1, inputfile[ict] );
 	}
 
@@ -1590,12 +1629,13 @@
 	
 	// read the direction of the incoming shower
-
+	// It is done only for one telescope since it is suposed
+	// to be the same shower for all of them
 	if (reflector_file_version<6){
-	  thetashw = mcevth.get_theta();
-	  phishw = mcevth.get_phi();
+	  thetashw = mcevth[0].get_theta();
+	  phishw = mcevth[0].get_phi();
 	}
 	else{
-	  thetashw = mcevth_2.get_theta();
-	  phishw = mcevth_2.get_phi();
+	  thetashw = mcevth_2[0].get_theta();
+	  phishw = mcevth_2[0].get_phi();
 	}
 	
@@ -1608,74 +1648,75 @@
 	if (reflector_file_version<6){
 
-	  mcevth.get_core(&coreX, &coreY);	
+	  mcevth[0].get_core(&coreX, &coreY);	
 	  
-	  // read the deviation of the telescope with respect to the shower
-	  mcevth.get_deviations ( &thetaCT, &phiCT );
-
-	  if ( (thetaCT == 0.) && (phiCT == 0.) ) {
+	  for(int ict=0;ict<ct_Number;ict++){
+	    // read the deviation of the telescope with respect to the shower
+	    mcevth[ict].get_deviations ( &thetaCT[ict], &phiCT[ict] );
+
+	    if ( (thetaCT[ict] == 0.) && (phiCT[ict] == 0.) ) {
 	  
-	    // CT was looking to the source (both lines are parallel)
-	    // therefore, we calculate the impact parameter as the distance 
-	    // between the CT axis and the core position
-
-	    impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
-	    thetaCT += thetashw;
-	    phiCT   += phishw;
+	      // CT was looking to the source (both lines are parallel)
+	      // therefore, we calculate the impact parameter as the distance 
+	      // between the CT axis and the core position
+
+	      impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
+	      thetaCT[ict] += thetashw;
+	      phiCT[ict]   += phishw;
+	    } else {
 	    
-	  } else {
+	      // the shower comes off-axis
 	    
-	    // the shower comes off-axis
+	      // obtain with this the final direction of the CT
+	      
+	      thetaCT[ict] += thetashw;
+	      phiCT[ict]   += phishw;
+	      
+	      // calculate vector for telescope
 	    
-	    // obtain with this the final direction of the CT
-	    
-	    thetaCT += thetashw;
-	    phiCT   += phishw;
-	    
-	    // calculate vector for telescope
-	    
-	    l2 = sin(thetaCT)*cos(phiCT);
-	    m2 = sin(thetaCT)*sin(phiCT);
-	    n2 = cos(thetaCT);
-	    
-	    num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
-	    den = (SQR(l1*m2 - l2*m1) + 
-		   SQR(m1*n2 - m2*n1) + 
-		   SQR(n1*l2 - n2*l1));
-	    den = sqrt(den);
-	    
-	    impactD = fabs(num)/den;
-	    
+	      l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
+	      m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
+	      n2 = cos(thetaCT[ict]);
+	      
+	      num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
+	      den = (SQR(l1*m2 - l2*m1) + 
+		     SQR(m1*n2 - m2*n1) + 
+		     SQR(n1*l2 - n2*l1));
+	      den = sqrt(den);
+	      impactD[ict] = fabs(num)/den;
+	    }
 	  }
 	}
 	else{
-	  mcevth_2.get_core(&coreX, &coreY);	
-	  thetaCT=mcevth_2.get_theta_CT();
-	  phiCT=mcevth_2.get_phi_CT();
-	  if ( (thetaCT == thetashw) && (phiCT == phishw) ) {
+	  mcevth_2[0].get_core(&coreX, &coreY);	
+	  for(int ict=0;ict<ct_Number;ict++){
+	    thetaCT[ict]=mcevth_2[ict].get_theta_CT();
+	    phiCT[ict]=mcevth_2[ict].get_phi_CT();
+	    if ( (thetaCT[ict] == thetashw) && (phiCT[ict] == phishw) ) {
 	  
-	    // CT was looking to the source (both lines are parallel)
-	    // therefore, we calculate the impact parameter as the distance 
-	    // between the CT axis and the core position
-
-	    impactD = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
+	      // CT was looking to the source (both lines are parallel)
+	      // therefore, we calculate the impact parameter as the distance 
+	      // between the CT axis and the core position
+
+	      impactD[ict] = dist_r_P( 0., 0., 0., l1, m1, n1, coreX, coreY, 0. );
 	    
-	  } else {
+	    } else {
 	    
-	    // the shower comes off-axis
+	      // the shower comes off-axis
 	    
-	    // calculate vector for telescope
+	      // calculate vector for telescope
 	    
-	    l2 = sin(thetaCT)*cos(phiCT);
-	    m2 = sin(thetaCT)*sin(phiCT);
-	    n2 = cos(thetaCT);
+	      l2 = sin(thetaCT[ict])*cos(phiCT[ict]);
+	      m2 = sin(thetaCT[ict])*sin(phiCT[ict]);
+	      n2 = cos(thetaCT[ict]);
 	    
-	    num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
-	    den = (SQR(l1*m2 - l2*m1) + 
-		   SQR(m1*n2 - m2*n1) + 
-		   SQR(n1*l2 - n2*l1));
-	    den = sqrt(den);
+	      num = (m1*n2*coreX - m2*n1*coreX + l2*n1*coreY - l1*n2*coreY);
+	      den = (SQR(l1*m2 - l2*m1) + 
+		     SQR(m1*n2 - m2*n1) + 
+		     SQR(n1*l2 - n2*l1));
+	      den = sqrt(den);
 	    
-	    impactD = fabs(num)/den;
+	      impactD[ict] = fabs(num)/den;
 	    
+	    }
 	  }
 	}
@@ -1685,12 +1726,12 @@
 	if ( Select_Energy ) {
 	  if (reflector_file_version<6)
-	    if (( mcevth.get_energy() < Select_Energy_le ) ||
-		( mcevth.get_energy() > Select_Energy_ue )) {
+	    if (( mcevth[0].get_energy() < Select_Energy_le ) ||
+		( mcevth[0].get_energy() > Select_Energy_ue )) {
 	      log(SIGNATURE, "select_energy: shower rejected.\n");
 	      continue;
 	    }
 	  else
-	    if (( mcevth_2.get_energy() < Select_Energy_le ) ||
-		( mcevth_2.get_energy() > Select_Energy_ue )) {
+	    if (( mcevth_2[0].get_energy() < Select_Energy_le ) ||
+		( mcevth_2[0].get_energy() > Select_Energy_ue )) {
 	      log(SIGNATURE, "select_energy: shower rejected.\n");
 	      continue;
@@ -1699,11 +1740,18 @@
 
 	// Read first and last time and put inumphe_CT[0] to 0
-
-	if (reflector_file_version<6)
-	  mcevth.get_times(&arrtmin_ns,&arrtmax_ns);
-	else
-	  mcevth_2.get_times(&arrtmin_ns,&arrtmax_ns);
-
-	ncph_system=0;
+	float tm=0, tM=9E+50;
+	for(int ict=0;ict<ct_Number;ict++){
+
+	  if (reflector_file_version<6)
+	    mcevth[ict].get_times(&arrtmin_ns,&arrtmax_ns);
+	  else
+	    mcevth_2[ict].get_times(&arrtmin_ns,&arrtmax_ns);
+	  
+	  tm=(tm<arrtmin_ns)?tm:arrtmin_ns;
+	  tM=(tM<arrtmax_ns)?tM:arrtmax_ns;
+	    
+	}
+	arrtmin_ns=tm;arrtmax_ns=tM;
+
 	inumphe=0;
 
@@ -1722,5 +1770,5 @@
 			    &inumphe_CT[ict], // important for later: the size of photoe[]
 			    fnpix,    // will be changed by the function!
-			    &ncph,    // will be changed by the function!
+			    &ncph[ict],    // will be changed by the function!
 			    &arrtmin_ns, // will be changed by the function!
 			    &arrtmax_ns, // will be changed by the function!
@@ -1728,5 +1776,5 @@
 
 	  inumphe=(inumphe<inumphe_CT[ict])?inumphe_CT[ict]:inumphe;
-	  
+
 	  if( k != 0 ){ // non-zero returnvalue means error
 	    cout << "Exiting.\n";
@@ -1795,22 +1843,26 @@
 	
 	}// end if(simulateNSB && nphe2NSB<=inumphe_CT[0]) ...
-  
-	inumphensb=0;
-	for (UInt_t ui=0;ui<
-	       ((MGeomCam*)(camgeom.UncheckedAt(0)))->GetNumPixels();ui++){
-	  inumphensb+=nsb_phepns[0][ui]*TOTAL_TRIGGER_TIME;
+
+	for(int ict=0;ict<ct_Number;ict++){
+	  inumphensb[ict]=0;
+	  for (UInt_t ui=0;ui<
+		 ((MGeomCam*)(camgeom.UncheckedAt(ict)))->GetNumPixels();ui++){
+	    inumphensb[ict]+=nsb_phepns[ict][ui]*TOTAL_TRIGGER_TIME;
+	  }
+	  ntcph[ict]+=ncph[ict];
+	  if ((nshow+ntshow)%100 == 1){
+	    log(SIGNATURE, "End of this event: %d cphs(+%d). . .\n",
+		ncph[ict], ntcph[ict]);
+
+	    cout << "Total number of phes in CT "<<ict<<": " 
+		 << inumphe_CT[ict]<<" (+ ";
+	    cout<<inumphensb[ict]<<" mean expected number from NSB)"<<endl;
+	  }
 	}
-	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 in CT 0: " << inumphe_CT[0]<<" + ";
-	  cout<<inumphensb<<endl;
-	}
 
 	// skip it ?
-	
-	for ( i=0; i<nSkip; ++i ) {
+
+	int i;	
+	for (i=0; i<nSkip; ++i ) {
 	  if (Skip[i] == (nshow+ntshow)) {
 	    i = -1;
@@ -1878,5 +1930,5 @@
 	  int iconcount;
 	  for (iconcount=0, ithrescount=0, fthrescount=Trigger_loop_lthres;fthrescount<=Trigger_loop_uthres;ithrescount++, fthrescount+=Trigger_loop_sthres){
-	    for (i=0;i<ct_NPixels;i++){
+	    for (int i=0;i<ct_NPixels;i++){
 	      fpixelthres[i]=
 		((Float_t)(fthrescount)>=qThreshold[0][i])?
@@ -1972,5 +2024,5 @@
 		      EvtHeader[iconcount]->FillHeader ( (UInt_t) (ntshow + nshow),0);
 		      //   fill pixel information
-		      if (Lev1){
+		      if (Lev1 || Write_All_Images){
 			if (addElecNoise) Fadc_CT[0]->DigitalNoise();
 			for(UInt_t i=0;
@@ -2013,19 +2065,19 @@
 	      Float_t ftime, ltime;
 	      if (reflector_file_version<6){
-		mcevth.get_times(&ftime, &ltime);
-		McEvt->Fill( 0, 
-			     (UShort_t) mcevth.get_primary() , 
-			     mcevth.get_energy(),
+		mcevth[0].get_times(&ftime, &ltime);
+		McEvt[0]->Fill( 0, 
+			     (UShort_t) mcevth[0].get_primary() , 
+			     mcevth[0].get_energy(),
 			     -1.0,
 			     -1.0,
 			     -1.0, 
-			     mcevth.get_theta(), 
-			     mcevth.get_phi(), 
-			     mcevth.get_core(),
-			     mcevth.get_coreX(),
-			     mcevth.get_coreY(),
-			     impactD,
-			     phiCT,
-			     thetaCT,
+			     mcevth[0].get_theta(), 
+			     mcevth[0].get_phi(), 
+			     mcevth[0].get_core(),
+			     mcevth[0].get_coreX(),
+			     mcevth[0].get_coreY(),
+			     impactD[0],
+			     phiCT[0],
+			     thetaCT[0],
 			     ftime,
 			     ltime,
@@ -2037,10 +2089,10 @@
 			     0,
 			     0,
-			     (UInt_t)mcevth.get_CORSIKA(), 
-			     (UInt_t)mcevth.get_AtmAbs(), 
-			     (UInt_t)(mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot()), 
-			     (UInt_t) ncph, 
+			     (UInt_t)mcevth[0].get_CORSIKA(), 
+			     (UInt_t)mcevth[0].get_AtmAbs(), 
+			     (UInt_t)(mcevth[0].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()), 
+			     (UInt_t) ncph[0], 
 			     (UInt_t) inumphe_CT[0],
-			     (UInt_t) inumphensb+inumphe_CT[0],
+			     (UInt_t) inumphensb[0]+inumphe_CT[0],
 			     -1.0,
 			     -1.0,
@@ -2049,20 +2101,20 @@
 	      else{
 		Float_t Nmax, t0, tmax, a, b, c, chi2;
-		mcevth_2.get_times(&ftime, &ltime);
-		chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
-		McEvt->Fill((UInt_t) mcevth_2.get_evt_number(),
-			    (UShort_t) mcevth_2.get_primary() , 
-			     mcevth_2.get_energy(),
-			     mcevth_2.get_thick0(),
-			     mcevth_2.get_first_target(),
-			     mcevth_2.get_z_first_int(),
-			     mcevth_2.get_theta(), 
-			     mcevth_2.get_phi(), 
-			     mcevth_2.get_core(),
-			     mcevth_2.get_coreX(),
-			     mcevth_2.get_coreY(),
-			     impactD,
-			     mcevth_2.get_phi_CT(),
-			     mcevth_2.get_theta_CT(),
+		mcevth_2[0].get_times(&ftime, &ltime);
+		chi2=mcevth_2[0].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
+		McEvt[0]->Fill((UInt_t) mcevth_2[0].get_evt_number(),
+			    (UShort_t) mcevth_2[0].get_primary() , 
+			     mcevth_2[0].get_energy(),
+			     mcevth_2[0].get_thick0(),
+			     mcevth_2[0].get_first_target(),
+			     mcevth_2[0].get_z_first_int(),
+			     mcevth_2[0].get_theta(), 
+			     mcevth_2[0].get_phi(), 
+			     mcevth_2[0].get_core(),
+			     mcevth_2[0].get_coreX(),
+			     mcevth_2[0].get_coreY(),
+			     impactD[0],
+			     mcevth_2[0].get_phi_CT(),
+			     mcevth_2[0].get_theta_CT(),
 			     ftime,
 			     ltime,
@@ -2074,13 +2126,13 @@
 			     c,
 			     chi2,
-			     (UInt_t)mcevth_2.get_CORSIKA(), 
-			     (UInt_t)mcevth_2.get_AtmAbs(), 
-			     (UInt_t)(mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot()), 
-			     (UInt_t) ncph, 
+			     (UInt_t)mcevth_2[0].get_CORSIKA(), 
+			     (UInt_t)mcevth_2[0].get_AtmAbs(), 
+			     (UInt_t)(mcevth_2[0].get_MirrAbs()+mcevth_2[0].get_OutOfMirr()+mcevth_2[0].get_BlackSpot()), 
+			     (UInt_t) ncph[0], 
 			     (UInt_t) inumphe_CT[0],
-			     (UInt_t) inumphensb+inumphe_CT[0],
-			     mcevth_2.get_ElecFraction(),
-			     mcevth_2.get_MuonFraction(),
-			     mcevth_2.get_OtherFraction());
+			     (UInt_t) inumphensb[0]+inumphe_CT[0],
+			     mcevth_2[0].get_ElecFraction(),
+			     mcevth_2[0].get_MuonFraction(),
+			     mcevth_2[0].get_OtherFraction());
 	      }
 	    }
@@ -2090,10 +2142,10 @@
 	    //  Clear the branches
 	    if(Write_McTrig){
-	      for(i=0;i<numberBranches;i++){
+	      for(int i=0;i<numberBranches;i++){
 		McTrig[i]->Clear() ;
 	      }
 	    }
 	    if( Write_RawEvt ){
-	      for(i=0;i<numberBranches;i++){
+	      for(int i=0;i<numberBranches;i++){
 		EvtHeader[i]->Clear() ;
 		EvtData[i]->DeletePixels();
@@ -2101,5 +2153,5 @@
 	    }
 	    if (Write_McEvt)
-	      McEvt->Clear() ;  
+	      McEvt[0]->Clear() ;  
 	  }
 	}
@@ -2140,5 +2192,5 @@
 	    
 	    Lev0MT[ict] = (Short_t) Trigger_CT[ict]->ZeroLevel() ; 
-	  
+
 	    Lev1MT[ict] = 0 ; 
 	  
@@ -2211,5 +2263,5 @@
 		//   fill pixel information
 	      
-		if (Lev1MT[ict]){
+		if (Lev1MT[ict] || Write_All_Images){
 		   if (addElecNoise) Fadc_CT[ict]->DigitalNoise();
 		  for(UInt_t i=0;
@@ -2234,5 +2286,5 @@
 	    if(FADC_Scan){
 	      if ( Lev0MT[ict] > 0 ) {
-		  Fadc_CT[ict]->ShowSignal( McEvt, (Float_t) 60. ) ; 
+		  Fadc_CT[ict]->ShowSignal( McEvt[ict], (Float_t) 60. ) ; 
 	      }
 	    }
@@ -2240,5 +2292,5 @@
 	    if(Trigger_Scan){
 	      if ( Lev0MT[ict] > 0 ) {
-		Trigger_CT[ict]->ShowSignal(McEvt) ;
+		Trigger_CT[ict]->ShowSignal(McEvt[ict]) ;
 	      }
 	    }
@@ -2248,81 +2300,83 @@
 	  // If there is trigger in some telecope or we store all showers
 	  if(btrigger){  
-	    //
-	    //   fill the MMcEvt with all information  
-	    //
+	    if (Write_McEvt){
+	      //
+	      //   fill the MMcEvt with all information  
+	      //
 	    
-	    if (Write_McEvt){
-	      Float_t ftime, ltime;
-	      if (reflector_file_version<6){
-		mcevth.get_times(&ftime, &ltime);
-		McEvt->Fill( 0, 
-			     (UShort_t) mcevth.get_primary() , 
-			     mcevth.get_energy(),
-			     -1.0,
-			     -1.0,
-			     -1.0, 
-			     mcevth.get_theta(), 
-			     mcevth.get_phi(), 
-			     mcevth.get_core(),
-			     mcevth.get_coreX(),
-			     mcevth.get_coreY(),
-			     impactD,
-			     phiCT,
-			     thetaCT,
-			     ftime,
-			     ltime,
-			     0,
-			     0,
-			     0,
-			     0,
-			     0,
-			     0,
-			     0,
-			     (UInt_t)mcevth.get_CORSIKA(), 
-			     (UInt_t)mcevth.get_AtmAbs(), 
-			     (UInt_t)(mcevth.get_MirrAbs()+mcevth.get_OutOfMirr()+mcevth.get_BlackSpot()), 
-			     (UInt_t) ncph, 
-			     (UInt_t) inumphe_CT[0],
-			     (UInt_t) inumphensb+inumphe_CT[0],
-			     -1.0,
-			     -1.0,
-			     -1.0);
-	      }
-	      else{
-		Float_t Nmax, t0, tmax, a, b, c, chi2;
-		mcevth_2.get_times(&ftime, &ltime);
-		chi2=mcevth_2.get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
-		McEvt->Fill( (UInt_t) mcevth_2.get_evt_number(),
-			     (UShort_t) mcevth_2.get_primary() , 
-			     mcevth_2.get_energy(),
-			     mcevth_2.get_thick0(),
-			     mcevth_2.get_first_target(),
-			     mcevth_2.get_z_first_int(),
-			     mcevth_2.get_theta(), 
-			     mcevth_2.get_phi(), 
-			     mcevth_2.get_core(),
-			     mcevth_2.get_coreX(),
-			     mcevth_2.get_coreY(),
-			     impactD,
-			     mcevth_2.get_phi_CT(),
-			     mcevth_2.get_theta_CT(),
-			     ftime,
-			     ltime,
-			     Nmax,
-			     t0,
-			     tmax,
-			     a,
-			     b,
-			     c,
-			     chi2,
-			     (UInt_t)mcevth_2.get_CORSIKA(), 
-			     (UInt_t)mcevth_2.get_AtmAbs(), 
-			     (UInt_t) (mcevth_2.get_MirrAbs()+mcevth_2.get_OutOfMirr()+mcevth_2.get_BlackSpot()), 
-			     (UInt_t) ncph, 
-			     (UInt_t) inumphe_CT[0],
-			     (UInt_t) inumphensb+inumphe_CT[0],
-			     mcevth_2.get_ElecFraction(),
-			     mcevth_2.get_MuonFraction(),
-			     mcevth_2.get_OtherFraction());
+	      for (int ict=0;ict<ct_Number;ict++){
+		Float_t ftime, ltime;
+		if (reflector_file_version<6){
+		  mcevth[ict].get_times(&ftime, &ltime);
+		  McEvt[ict]->Fill( 0, 
+				  (UShort_t) mcevth[0].get_primary() , 
+				  mcevth[ict].get_energy(),
+				  -1.0,
+				  -1.0,
+				  -1.0, 
+				  mcevth[ict].get_theta(), 
+				  mcevth[ict].get_phi(), 
+				  mcevth[ict].get_core(),
+				  mcevth[ict].get_coreX(),
+				  mcevth[ict].get_coreY(),
+				  impactD[ict],
+				  phiCT[ict],
+				  thetaCT[ict],
+				  ftime,
+				  ltime,
+				  0,
+				  0,
+				  0,
+				  0,
+				  0,
+				  0,
+				  0,
+				  (UInt_t)mcevth[ict].get_CORSIKA(), 
+				  (UInt_t)mcevth[ict].get_AtmAbs(), 
+				  (UInt_t)(mcevth[ict].get_MirrAbs()+mcevth[0].get_OutOfMirr()+mcevth[0].get_BlackSpot()), 
+				  (UInt_t) ncph[ict], 
+				  (UInt_t) inumphe_CT[ict],
+				  (UInt_t) inumphensb[ict]+inumphe_CT[ict],
+				  -1.0,
+				  -1.0,
+				  -1.0);
+		}
+		else{
+		  Float_t Nmax, t0, tmax, a, b, c, chi2;
+		  mcevth_2[ict].get_times(&ftime, &ltime);
+		  chi2=mcevth_2[ict].get_NKGfit(&Nmax, &t0, &tmax, &a, &b, &c);
+		  McEvt[ict]->Fill( (UInt_t) mcevth_2[ict].get_evt_number(),
+				  (UShort_t) mcevth_2[ict].get_primary() , 
+				  mcevth_2[ict].get_energy(),
+				  mcevth_2[ict].get_thick0(),
+				  mcevth_2[ict].get_first_target(),
+				  mcevth_2[ict].get_z_first_int(),
+				  mcevth_2[ict].get_theta(), 
+				  mcevth_2[ict].get_phi(), 
+				  mcevth_2[ict].get_core(),
+				  mcevth_2[ict].get_coreX(),
+				  mcevth_2[ict].get_coreY(),
+				  impactD[ict],
+				  mcevth_2[ict].get_phi_CT(),
+				  mcevth_2[ict].get_theta_CT(),
+				  ftime,
+				  ltime,
+				  Nmax,
+				  t0,
+				  tmax,
+				  a,
+				  b,
+				  c,
+				  chi2,
+				  (UInt_t)mcevth_2[ict].get_CORSIKA(), 
+				  (UInt_t)mcevth_2[ict].get_AtmAbs(), 
+				  (UInt_t) (mcevth_2[ict].get_MirrAbs()+mcevth_2[ict].get_OutOfMirr()+mcevth_2[ict].get_BlackSpot()), 
+				  (UInt_t) ncph[ict], 
+				  (UInt_t) inumphe_CT[ict],
+				  (UInt_t) inumphensb[ict]+inumphe_CT[ict],
+				  mcevth_2[ict].get_ElecFraction(),
+				  mcevth_2[ict].get_MuonFraction(),
+				  mcevth_2[ict].get_OtherFraction());
+		}
 	      }
 	    }
@@ -2342,6 +2396,6 @@
 	    if (Write_RawEvt) EvtData[ict]->DeletePixels();
 	    if (Write_McTrig) McTrig[ict]->Clear() ;
+	    if (Write_McEvt) McEvt[ict]->Clear() ; 
 	  } 
-	  if (Write_McEvt) McEvt->Clear() ; 
 	}
 		
@@ -2380,5 +2434,5 @@
 
 	// We search the maximum impact parameter fo the simualted showers
-	maxpimpact=maxpimpact<impactD?impactD:maxpimpact;
+	maxpimpact=maxpimpact<impactD[0]?impactD[0]:maxpimpact;
 	
 	// look for the next event
@@ -2442,23 +2496,23 @@
     telestheta=-10.0;
     telesphi=-10.0;
-    mcevth.get_theta_range(&shthetamin, &shthetamax);
-    mcevth.get_phi_range(&shphimin,&shphimax);
-    mcevth.get_theta_range(&shthetamin, &shthetamax);
-    mcevth.get_phi_range(&shphimin,&shphimax);
-    corsika=UInt_t(mcevth.get_VersionPGM()*1000);
+    mcevth[0].get_theta_range(&shthetamin, &shthetamax);
+    mcevth[0].get_phi_range(&shphimin,&shphimax);
+    mcevth[0].get_theta_range(&shthetamin, &shthetamax);
+    mcevth[0].get_phi_range(&shphimin,&shphimax);
+    corsika=UInt_t(mcevth[0].get_VersionPGM()*1000);
     for (int i=0; i< 10;i++)
-      heights[i]=mcevth.get_HeightLev (i);
-    rnum=mcevth.get_RunNumber();
+      heights[i]=mcevth[0].get_HeightLev (i);
+    rnum=mcevth[0].get_RunNumber();
   }
   else{
-    telestheta=mcevth_2.get_theta_CT();
-    telesphi=mcevth_2.get_phi_CT();
-    mcevth_2.get_theta_range(&shthetamin, &shthetamax);
-    mcevth_2.get_phi_range(&shphimin,&shphimax);
-    corsika=UInt_t(mcevth_2.get_VersionPGM()*1000);
+    telestheta=mcevth_2[0].get_theta_CT();
+    telesphi=mcevth_2[0].get_phi_CT();
+    mcevth_2[0].get_theta_range(&shthetamin, &shthetamax);
+    mcevth_2[0].get_phi_range(&shphimin,&shphimax);
+    corsika=UInt_t(mcevth_2[0].get_VersionPGM()*1000);
     for (int i=0; i< 10;i++)
-      heights[i]=mcevth_2.get_HeightLev (i);
-    rnum=mcevth_2.get_RunNumber();
-    mcevth_2.get_viewcone(&viewcone[0],&viewcone[1]);
+      heights[i]=mcevth_2[0].get_HeightLev (i);
+    rnum=mcevth_2[0].get_RunNumber();
+    mcevth_2[0].get_viewcone(&viewcone[0],&viewcone[1]);
   }
 
@@ -2471,5 +2525,5 @@
     McRunHeader->Fill(rnum,
 		    (UInt_t) 0,
-		    mcevth.get_DateRun(),
+		    mcevth[0].get_DateRun(),
 		    ftime,
 		    icontrigger,
@@ -2500,7 +2554,7 @@
 		    shphimin,
 		    maxpimpact,
-		    mcevth.get_CWaveLower(),
-		    mcevth.get_CWaveUpper(),
-		    mcevth.get_slope(),
+		    mcevth[0].get_CWaveLower(),
+		    mcevth[0].get_CWaveUpper(),
+		    mcevth[0].get_slope(),
 		    1,
 		    heights,
@@ -2512,5 +2566,5 @@
     McRunHeader->Fill(rnum,
 		    (UInt_t) 0,
-		    mcevth_2.get_DateRun(),
+		    mcevth_2[0].get_DateRun(),
 		    ftime,
 		    icontrigger,
@@ -2541,7 +2595,7 @@
 		    shphimin,
 		    maxpimpact,
-		    mcevth_2.get_CWaveLower(),
-		    mcevth_2.get_CWaveUpper(),
-		    mcevth_2.get_slope(),
+		    mcevth_2[0].get_CWaveLower(),
+		    mcevth_2[0].get_CWaveUpper(),
+		    mcevth_2[0].get_slope(),
 		    1,
 		    heights,
@@ -2580,5 +2634,5 @@
 			     1,
 			     heights,
-			     mcevth_2.get_slope(),
+			     mcevth_2[0].get_slope(),
 			     mcrunh.get_ELow(),
 			     mcrunh.get_EUpp(),
@@ -2638,8 +2692,8 @@
   // close input file
   
-  log( SIGNATURE, "%d event(s), with a total of %d C.photons\n", 
-       ntshow, ntcph );
-  datafile<<ntshow<<" event(s), with a total of "<<ntcph<<" C.photons"<<endl;
   if (Trigger_Loop){
+    log( SIGNATURE, "%d event(s), with a total of %d C.photons\n", 
+	 ntshow, ntcph[0] );
+    datafile<<ntshow<<" event(s), with a total of "<<ntcph[0]<<" C.photons"<<endl;
     log( SIGNATURE, "Trigger Mode. \n");
     log( SIGNATURE, "Fraction of triggers: \n");
@@ -2658,4 +2712,9 @@
   else{
     for(int ict=0;ict<ct_Number;ict++){
+      log( SIGNATURE, 
+	   "%d event(s), with a total of %d C.photons in CT %i (%s)\n", 
+	   ntshow, ntcph[ict],ict,GeometryName[ict] );
+      datafile<<ntshow<<" event(s), with a total of "<<ntcph[ict]
+	      <<" C.photons in CT "<<ict<<" ("<<GeometryName[ict]<<")"<<endl;
       log( SIGNATURE, "Fraction of triggers: %5.1f%% (%d out of %d)\n", 
 	   ((float)ntrigger[ict]) / ((float)ntshow) * 100.0, ntrigger[ict], ntshow);
@@ -3839,7 +3898,4 @@
   int k, ii; // counters
 
-  MTrigger trigger(camgeom->GetNumPixels(),camgeom,Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
-  MFadc flashadc(camgeom->GetNumPixels());
-
   static float wl_nm[iNUMWAVEBANDS + 1] = { WAVEBANDBOUND1,
 					    WAVEBANDBOUND2,
@@ -3886,4 +3942,9 @@
     exit(1);
   }
+
+  // Instance of MTrigger and MFadc
+
+  MTrigger trigger(camgeom->GetNumPixels(),camgeom,Trigger_gate_length, Trigger_overlaping_time, Trigger_response_ampl, Trigger_response_fwhm);
+  MFadc flashadc(camgeom->GetNumPixels());
 
   // initialize flag
@@ -4157,4 +4218,8 @@
 //
 // $Log: not supported by cvs2svn $
+// Revision 1.64  2003/10/21 07:42:50  blanch
+// A factor 2.35 to transform the fwhm into the sigma of gaussian was missing
+// in the storing of FADC single hpe pulse determination.
+//
 // Revision 1.63  2003/10/17 19:38:31  blanch
 // Now the camera program will stop if a undefined Geometry is required.