Changeset 2427 for trunk/MagicSoft/Simulation/Detector

Timestamp:

10/26/03 19:43:00 (21 years ago)

Author:

blanch

Message:

- The screen output information has been improved to prevent some
non-desired running conditions just looking at the output screen.
- One MMcEvt for each Telscope is stored in the output file.
- 500 empty events are simualted to get a more precise estimation of the
pedestal Sigma for each pixel.

File:

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

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

@@ -21 +21 @@
 //
 // $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 +165 @@
 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 +339 @@
 
   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 +361 @@
   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 +387 @@
   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 +453 @@
 
   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 +464 @@
       nsbrate_phepns[i][j]=0.0;    //@< Starlight rates initialised at 0.0
   }
+  for(int i=0;i<100;i++)
+    ntcph[i]=0;
   /*!@'
 
@@ -529 +528 @@
   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 +602 @@
   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 +747 @@
   
   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 +810 @@
 
     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 +885 @@
     McTrig = new MMcTrig * [numberBranches];
   
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       McTrig[i] = new MMcTrig();
     }
@@ -896 +891 @@
     HeaderTrig = new MMcTrigHeader * [numberBranches];
   
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       HeaderTrig[i] = new MMcTrigHeader();
     }
@@ -905 +900 @@
     HeaderFadc = new MMcFadcHeader * [numberBranches];
   
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       HeaderFadc[i] = new MMcFadcHeader();
     }
@@ -934 +929 @@
   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 +956 @@
   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 +967 @@
     EvtHeader = new MRawEvtHeader * [numberBranches]; 
 
-    for (i=0;i<numberBranches;i++) {
+    for (int i=0;i<numberBranches;i++) {
       EvtHeader[i] = new MRawEvtHeader();
     }
@@ -984 +979 @@
     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 +986 @@
   }
 
-  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 +1054 @@
 
   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 +1072 @@
   }
   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 +1147 @@
   //  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 +1167 @@
     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 +1183 @@
       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 +1204 @@
           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 +1223 @@
   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 +1258 @@
 
     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 +1272 @@
 
   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 +1334 @@
     
     // FIXME --- star NSB different for each camera?
+    log(SIGNATURE,"Produce NSB rates from Star Field");
+
     k = produce_nsbrates( starfieldname,
                           ((MGeomCam*)(camgeom.UncheckedAt(0))),
@@ -1352 +1381 @@
         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 +1404 @@
   //!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   //
-  //  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 +1414 @@
 
   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 +1591 @@
         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 +1629 @@
         
         // 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 +1648 @@
         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 +1726 @@
         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 +1740 @@
 
         // 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 +1770 @@
                             &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 +1776 @@
 
           inumphe=(inumphe<inumphe_CT[ict])?inumphe_CT[ict]:inumphe;
-          
+
           if( k != 0 ){ // non-zero returnvalue means error
             cout << "Exiting.\n";
@@ -1795 +1843 @@
         
         }// 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 +1930 @@
           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 +2024 @@
                       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 +2065 @@
               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 +2089 @@
                              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 +2101 @@
               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 +2126 @@
                              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 +2142 @@
             //  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 +2153 @@
             }
             if (Write_McEvt)
-              McEvt->Clear() ;  
+              McEvt[0]->Clear() ;  
           }
         }
@@ -2140 +2192 @@
             
             Lev0MT[ict] = (Short_t) Trigger_CT[ict]->ZeroLevel() ; 
-          
+
             Lev1MT[ict] = 0 ; 
           
@@ -2211 +2263 @@
                 //   fill pixel information
               
-                if (Lev1MT[ict]){
+                if (Lev1MT[ict] || Write_All_Images){
                    if (addElecNoise) Fadc_CT[ict]->DigitalNoise();
                   for(UInt_t i=0;
@@ -2234 +2286 @@
             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 +2292 @@
             if(Trigger_Scan){
               if ( Lev0MT[ict] > 0 ) {
-                Trigger_CT[ict]->ShowSignal(McEvt) ;
+                Trigger_CT[ict]->ShowSignal(McEvt[ict]) ;
               }
             }
@@ -2248 +2300 @@
           // 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 +2396 @@
             if (Write_RawEvt) EvtData[ict]->DeletePixels();
             if (Write_McTrig) McTrig[ict]->Clear() ;
+            if (Write_McEvt) McEvt[ict]->Clear() ; 
           } 
-          if (Write_McEvt) McEvt->Clear() ; 
         }
                 
@@ -2380 +2434 @@
 
         // 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 +2496 @@
     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 +2525 @@
     McRunHeader->Fill(rnum,
                     (UInt_t) 0,
-                    mcevth.get_DateRun(),
+                    mcevth[0].get_DateRun(),
                     ftime,
                     icontrigger,
@@ -2500 +2554 @@
                     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 +2566 @@
     McRunHeader->Fill(rnum,
                     (UInt_t) 0,
-                    mcevth_2.get_DateRun(),
+                    mcevth_2[0].get_DateRun(),
                     ftime,
                     icontrigger,
@@ -2541 +2595 @@
                     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 +2634 @@
                              1,
                              heights,
-                             mcevth_2.get_slope(),
+                             mcevth_2[0].get_slope(),
                              mcrunh.get_ELow(),
                              mcrunh.get_EUpp(),
@@ -2638 +2692 @@
   // 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 +2712 @@
   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 +3898 @@
   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 +3942 @@
     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 +4218 @@
 //
 // $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.