//////////////////////////////////////////////////////////////////
//
// checkmc                  
//
// @file        checkmc.cxx
// @title       Check of the MC generated data
// @author      J C Gonz\'alez, P. Morawitz, D. Petry
// @email       gonzalez@mppmu.mpg.de, morawitz@cern.ch, petry@ifae.es
//
// @maintitle
//_______________________________________________________________
//
//  Created: Thu May  7 16:24:22 1998
//  Author:  Jose Carlos Gonzales et al.
//  Purpose: Check of the MC generated data
//  Notes:  version 27-10-98 
//    
//    

// Below you can find a sort of, incomplete for sure, documentation 
// concerning this program. Please, do not hesiate to contact the 
// author in case of problems.

// @T \newpage

// @section Source code of {\tt checkmc.cxx}

/* @text
In this section we show the (commented) code of the program 
for for the check of the MC generated data, in the current version.
@endtext */

// @subsection Includes and Global variables definition

/* @text 
All the defines are located in the file {\tt checkmc.h}.
@endtext */
 
// @code
#include "checkmc.h"
// @endcode

// @subsubsection Definition of global variables.

/* @text
Here we define the global variables where the values from
the parameters file are stored. See the section
\ref{paramdesc} in page \pageref{paramdesc} for a more 
complete description of the available commands in the 
parameters file.
@endtext */

// @code
static char **Paths_list;     // list of paths to get data from
static int Num_of_paths = 1;  // number of paths
static char Output_filename[PATH_MAX_LENGTH]; // output filename
static char CT_filename[PATH_MAX_LENGTH]; // name of the CT def. file
static VerboseLevel verbose = VERBOSE_DEFAULT; // verboseness
static float fixed_Theta;      // fixed target theta and phi
static float fixed_Phi;
static float low_Ecut = 0.0, high_Ecut = 100000.0;
static int is_Fixed_Target = FALSE; // are we using fixed target?
static int max_Events = 100000; // maximum number of events to read

static float gammas  =    1.;
static float protons =   14.;
static float helium =   402.;


inline int max(int a, int b)
{if (a>b) return a;
 else return b;}

inline int min(int a, int b)
{if (a<b) return a;
 else return b;}

// @endcode

/* @text
Earth's radius (in cm).
@endtext */

// @code
static const float REarth = 6348.0e5;
// @endcode

// @subsection Main program

/* @text
Let's start!
@endtext */

// @code
//++++++++++++++++++++++++++++++++++++++++
// MAIN PROGRAM 
//----------------------------------------
                   
void
main(int argc, char **argv) 
{

  FILE *f;
  DIR *directory;             // directory of data
  struct dirent *de;          // directory entry

  char pathname[256];         // directory and file names
  char cername[256];
  char staname[256];
  char outname[256];
  char refname[256];
  char name[256];

  ifstream cerfile;
  ofstream outputfile;

  COREventHeader evth;        // Event Header class (from CORSIKA)
  CORParticle photon;         // Particle (photon) class (from CORSIKA)
  CORStatfile stat;           // Statistics file class (from CORSIKA)

  float thetaCT, phiCT, xiCT; // parameters of a given shower
  float coreD, coreX, coreY;  // core position and distance

  float a, b, c, t, tt;       // intermediate variables

  float mind, d;              // minimum distance trayectory-mirror center

  float wl;                   // wavelength of the Cphoton
  float reflec;               // reflectivity for a Cphoton

  int ncer;                   // number of the shower;
  int np, nf;                 // number of path, number of cer-file
  int i, j, k, n;             // simple counters
  int i_mirror=-1;            // number of a given mirror

  int nCphotons;              // number of a Cphotons written on output file

  float TR;                   // atmospheric transmitance
  int nbeforeTR, nafterTR;    // number of Cph. before and after transmission
  int num_cer_files;          // number of cer files;
  int max_num_cer_files;      // maximum number of cer files to read

  float t1, t2;               // first and last time

  float x0, xa, xb, v;        // variables for the 
  float vmax, vmaxa, vmaxb;   // estimation of Xmax
  float xmax;

  int ntotfiles;
  int kbin;

  int  hprimary;
  static const int nprimaries=4;
  char* primary_name[nprimaries] ={ 
       "Photons",
       "Protons",
       "Helium",
       "Heavy Nuclei"};
  int nevents[nprimaries];
  char htit[256]; 

  int firstevent,lastevent;
  int ndiag=-1;
  int rndiag=-1;
  char diag[256][20];
  char rdiag[256][100];
  char ctmp[256];       
  float par[3],vp2[4],vp3[5],vp4[6];
  float sigpar;
  float chi2;
  float rnd[3][3];

  // now i/p the constants to be checked
  const float mon_HeightLev=220000.;
  const float mon_SlopeSpec[nprimaries] =  {-1.5,-2.75,-2.62,-100000.};
  const float mon_ElowLim [nprimaries] = {1.,10.,10.,10.};
  const float mon_EUppLim[nprimaries] = {30000.,100000.,100000.,100000.};
  const float mon_Ecutoffh[nprimaries] = {0.3,0.3,0.3,0.3};
  const float mon_Ecutoffm[nprimaries] = {0.3,0.3,0.3,0.3};
  const float mon_Ecutoffe[nprimaries] = {0.02,0.02,0.02,0.02};
  const float mon_Ecutoffg[nprimaries] = {0.02,0.02,0.02,0.02};
  const float mon_EGS4yn[nprimaries] = {1.,1.,1.,1.};
//  the EGS4yn flag should be automatically 
  const float mon_NKGyn[nprimaries] = {0.,0.,0.,0.};
  const float mon_GHEISHAyn[nprimaries] = {1.,1.,1.,1.};
  const float mon_VENUSyn[nprimaries] = {1.,1.,1.,1.};
  const float mon_CERENKOVyn[nprimaries] = {1.,1.,1.,1.};
  const float mon_ThetaMin[nprimaries] = {5.,0.,0.,0.};
  const float mon_ThetaMax[nprimaries] = {25.,30.,30.,30.};
  const float mon_PhiMin[nprimaries] = {0.,0.,0.,0.};
  const float mon_PhiMax[nprimaries] = {360.,360.,360.,360.};
  const float mon_CWaveLower[nprimaries] = {290.,290.,290.,290.};
  const float mon_CWaveUpper[nprimaries] = {600.,600.,600.,600.};

  // @endcode

  // @text
  // Variables for the primary particle N-Tuple 
  // @endtext

  // @code
  int hid=1,istat=0,icycle=0;
  const int nvar=44;
  char chTags[nvar][8]={
	"primary",
//  primary particle id 
	"energy",
//  primary energy
	"cored",
//  Core distance  in cm
	"theta",
	"phi",
	"thick",
//  thickness (g/cm-2), e.g. height z converted, see below
	"height",
//  height, z coordinate, cm, of first interaction
	"target1",
//  number of first taget ?
	"px",
	"py",
	"pz",
//  the momentum of the primary particle
	"tfirst", 
//  time of first cerenkov photon on ground
	"tlast",   
//  time of last cerenkov photon on ground
	"xmax", 
//  xmax in g/cm**2; "height" of cerenkov shower maximum
//  => xmax/dichte == penetration depth
//  depth defined as == 120km=0 
// (120km - penetration depth - observation level (2km)) == real height of interaction
// 
	"xcore",
//  x position of core on ground in cm
	"ycore",
//  y position of core on ground in cm
        "Nph", 
//  number of cerenkov photons 
        "RunNum",
        "EvtNum",
        "DateRun",
        "Version",
        "HeigLev",
        "ESlope",
        "ELowLim",
        "EUppLim",
//
        "Ecutofh",
        "Ecutofm",
        "Ecutofe",
        "Ecutofg",
//      
        "EGS4yn",
        "NKGyn",
        "GHEISHA",
        "VENUS",
        "CERENKO",
        "NEUTRIN",
        "HORIZON",
        "COMPUTE",
        "ThetMin",
        "ThetMax",
        "PhiMin",
        "PhiMax",
//  
        "StepLen",
        "CWavLow",
        "CWavUpp"
  };
  // char chTags[nvar * 10];
  float rec[nvar];
  int record_size=4096;
  // @endcode


  // @text
  // Variables for the cerenkov photon N-Tuple 
  // @endtext

  // @code
  int hid2=2,istat2=0,icycle2=0;
  const int nvar2=11;
  char chTags2[nvar2][8]={
//
// First info about primary
	"energy",
	"cored",
	"theta",
//
// Now cerenkov photon info
	"Wavel", 
// wavelength in nm
        "xph",
        "yph",
        "uph",
        "vph",
        "tph",
        "hph",
        "EvtNum"
  };
  // char chTags2[nvar2 * 10];
  float rec2[nvar2];
  int record_size2=4096;
  // @endcode


  // @text
  // We start in the main program. First we (could) make some presentation, 
  // and follows the reading of the parameters file (now from the {\tt stdin}), 
  // the reading of the CT parameters file, and the creation of the 
  // output file, where the processed data will be stored.
  // @endtext
  // @code

  //++
  //  START
  //--

  // make unbuffered output

  cout.setf ( ios::stdio );

  // HBOOK core memory

  HLIMIT(PAWC_SIZE);

  /*
	strcpy( chTags, "primary:" );
	strcat( chTags, "energy:" );
	strcat( chTags, "cored:" );
	strcat( chTags, "theta:" );
	strcat( chTags, "phi:" );
	strcat( chTags, "thick:" );
	strcat( chTags, "height:" );
	strcat( chTags, "target1:" );
	strcat( chTags, "px:" );
	strcat( chTags, "py:" );
	strcat( chTags, "pz:" );
	strcat( chTags, "tfirst:" );
	strcat( chTags, "tlast:" );
	strcat( chTags, "xmax" );
  */

  // make some sort of presentation

  present();
  
  // read parameters from the stdin

  readparam();
  if ( verbose >= VERBOSE_MINIMAL )
     log( SIGNATURE, "... reading parameters OK\n");
  // get name of the output file
  
  strcpy( outname, get_output_filename() );

  // open HBOOK file

  if ( verbose >= VERBOSE_MINIMAL ) 
    log( SIGNATURE, "Openning N-Tuple file %s\n", outname );

  // TFile hfile(outname,"RECREATE","Check MC generated data");


  HROPEN(1, "CHECKMC", outname, "N", record_size, istat);
  if (istat != 0) 
      {fprintf(stderr, "\n Output Histo couldn't be opened. Adios.\n");
       exit(999);}

  // profile (final) histograms

  // TH1F *htime  = new TH1F("htime", "htime", 60, 0., 60.);
  // TH1F *helec  = new TH1F("helec", "helec", NTHSTEP, 0., 790.);
  // TH1F *hlight = new TH1F("hlight", "hlight", NTHSTEP, 0., 790.);

  double fstime[18][60], fs2time[18][60];

  double ftstime[18][60], fts2time[18][60];
  double ftselec[18][NTHSTEP], fts2elec[18][60];
  double ftslight[18][NTHSTEP], fts2light[18][60];
  
  for (i=0; i<60; ++i)
    for (k=0; k<18; ++k)
      ftstime[k][i] = fts2time[k][i] = 0.0;
  
  for (i=0; i<NTHSTEP; ++i)
    for (k=0; k<18; ++k)
      ftselec[k][i] = fts2elec[k][i] = 
		ftslight[k][i] = fts2light[k][i] = 0.0;
  
  // @endcode 
  // @text
  // After reading the parameters file (from {\tt stdin}, and the 
  // CT definition file, we start the analysis. Basically, there are  
  // three loops. The outermost loop is a loop in the {\it data directories\/} 
  // we introduced in the parameters file; the middle loop follows each 
  // {\tt cer*} file in the directory; and the innermost loop gets 
  // each Cherenkov from these files.
  // @endtext 
  // @code

  // TNtuple *ntuple = new TNtuple("ntuple","Check MC data", chTags);
  // HBOOKN(hid,  pathname, nvar,  "CHECKMC", 4096, chTags);
  // HBOOKN(hid2, pathname, nvar2, "CHECKMC", 4096, chTags2);



  // First count the number of events for each primary 
  ntotfiles = 0;
  // for each path of data files
  for (np=0; np<get_num_of_paths(); np++) {
    strcpy(pathname, get_path_name(np));
    // open directory 
    if ( verbose >= VERBOSE_MINIMAL )
      log( SIGNATURE, "Openning directory %s\n", pathname );
    if ( (directory = opendir(pathname)) == NULL ) 
      error( SIGNATURE, 
             "Cannot open directory %s\n", pathname );
    // trace the number of events (cer* files) 
    num_cer_files = 0;
    max_num_cer_files = get_max_events();
    // for each directory entry (files)
    while ( ( (de = readdir( directory )) != NULL) &&
            (num_cer_files < max_num_cer_files) ) {
      // skip removed files
      if ( de->d_ino == 0)
        continue;
      // keep only cer* files
      if ( strstr(de->d_name, "cer") != de->d_name )
        continue;
      // it is a real cer* file
      ++num_cer_files;
      if ( verbose >= VERBOSE_NORMAL )
        log( SIGNATURE, "cerfile: %s\n", de->d_name );
      // get cer* file number (number of the shower)
      // increments the pointer in 3 to calculate the number
      ncer = atoi(de->d_name + 3);  
      // full names of the input files cer* and sta*
      sprintf(cername, "%s/cer%06d", pathname, ncer);
      sprintf(staname, "%s/sta%06d", pathname, ncer);
      // try to open cer* file
      if ( verbose >= VERBOSE_MAXIMAL )
        log( SIGNATURE, "Openning %s\n", cername );
      cerfile.open( cername );
      if ( cerfile.bad() ) 
        error( SIGNATURE, 
               "Cannot open input file: %s\n", cername );
      // try to open the statistics file
      if ( verbose >= VERBOSE_MAXIMAL )
        log( SIGNATURE, "Openning %s\n", staname );
      stat.openfile ( staname );
      // reading data from this sta* file
      if ( verbose >= VERBOSE_MAXIMAL )
        log( SIGNATURE, "Reading %s\n", staname );
      stat.read();
      // read the Event Header from the cer* file
      evth.read( cerfile );
      if (evth.PrimaryID == gammas)  
         ++nevents[0];
      else if (evth.PrimaryID == protons)  
         ++nevents[1];
      else if (evth.PrimaryID == helium)  
         ++nevents[2];
      else
         ++nevents[3];
      cerfile.close();
      stat.closefile();}}


  cout<<"NMAXBINS"<<NMAXBINS<<endl;
  // Now book the histogram files  
  for (i=0; i<nprimaries; ++i)
      if (nevents[i]>0) 
// Primary histos
      {sprintf(htit, "Log(Energy) of primary %s;log(E/GeV)", primary_name[i]);
       HBOOK1(100+i,htit,min(nevents[i]/25,NMAXBINS),0.,4.5,0.); 

       sprintf(htit, "Core distance (horizontal plane) of primary %s;(m)", primary_name[i]);
       HBOOK1(110+i,htit,min(nevents[i]/25,NMAXBINS),0.,500.,0.); 

       sprintf(htit, "Theta of primary %s;(degrees)", primary_name[i]);
       HBOOK1(120+i,htit,min(nevents[i]/15,NMAXBINS),0.,70.,0.); 

       sprintf(htit, "Phi of primary %s;(degrees)", primary_name[i]);
       HBOOK1(130+i,htit,min(nevents[i]/25,NMAXBINS),0.,360.,0.); 

       sprintf(htit, "Core position for primary %s;x (m);y (m)", primary_name[i]);
       HBOOK2(140+i,htit,50,-500.,500.,50,-500.,500.,0.); 

       sprintf(htit, "Height of first interaction of primary %s;(km)", primary_name[i]);
       HBOOK1(150+i,htit,min(nevents[i]/25,NMAXBINS),0.,100.,0.); 
       
       sprintf(htit, "Time (last-first) of C-Photons for %s;(ns)", primary_name[i]);
       HBOOK1(160+i,htit,min(nevents[i]/25,NMAXBINS),-1.,160.,0.); 

       sprintf(htit, "Shower maximum XMAX for %s;(g cm^-2!)", primary_name[i]);
       HBOOK1(170+i,htit,min(nevents[i]/10,NTHSTEP),0.,790.,0.); 

       sprintf(htit, "Log(Num of C-photons produced) for primary %s below 300 GeV;log(N)", primary_name[i]);
       HBOOK1(180+i,htit,min(nevents[i]/25,NMAXBINS),-.5,6.,0.); 

       sprintf(htit, "Log(Num of C-photons produced) for primary %s above 300GeV;log(N)", primary_name[i]);
       HBOOK1(190+i,htit,min(nevents[i]/25,NMAXBINS),-.5,6.,0.); 

// C-Photon histos
       sprintf(htit, "C-Photon position in horiz. plane for prim. %s;x (m);y (m)", primary_name[i]);
       HBOOK2(300+i,htit,50,-15.,15.,50,-15.,15.,0.); 

       sprintf(htit, "C-Photon distance to (0,0) (horizontal plane) for %s;(m)", primary_name[i]);
       HBOOK1(310+i,htit,(min(nevents[i]/5,NMAXBINS)),0.,15.,0.); 

       sprintf(htit, "C-Photon spectrum for primary %s; [l] (nm)", primary_name[i]);
       HBOOK1(320+i,htit,(min(nevents[i]/5,NMAXBINS)),200.,700.,0.); 

       sprintf(htit, "Production height of C-Photon for primary %s;(km)", primary_name[i]);
       HBOOK1(330+i,htit,(min(nevents[i]/5,NMAXBINS)),0.,30.,0.); 
      }

  ntotfiles = 0;

  // for each path of data files
  
  for (np=0; np<get_num_of_paths(); np++) {
    
    strcpy(pathname, get_path_name(np));
    
    // open directory 

    if ( verbose >= VERBOSE_MINIMAL )
      log( SIGNATURE, "Openning directory %s\n", pathname );
    if ( (directory = opendir(pathname)) == NULL ) 
      error( SIGNATURE, 
             "Cannot open directory %s\n", pathname );
    
    // trace the number of events (cer* files) 
    
    num_cer_files = 0;

    // book n-tuple for this directory

    //  hid = np+1;

    //  if ( verbose >= VERBOSE_MINIMAL ) 
    //    log( SIGNATURE, "Create N-Tuple #%d for path %s . . .\n", 
    //             hid, pathname );
        
    // HBOOKN(hid, pathname, nvar, " ", 4096, chTags);
        
    max_num_cer_files = get_max_events();

    // for each directory entry (files)

    firstevent=-100;    
    while ( ( (de = readdir( directory )) != NULL) &&
            (num_cer_files < max_num_cer_files) ) {
	  
      // skip removed files
          
      if ( de->d_ino == 0)
        continue;

      // keep only cer* files

      if ( strstr(de->d_name, "cer") != de->d_name )
        continue;

      // it is a real cer* file

      ++num_cer_files;

      if ( verbose >= VERBOSE_NORMAL )
        log( SIGNATURE, "cerfile: %s\n", de->d_name );

      // get cer* file number (number of the shower)

      // increments the pointer in 3 to calculate the number

      ncer = atoi(de->d_name + 3);  
      
      // full names of the input files cer* and sta*

      sprintf(cername, "%s/cer%06d", pathname, ncer);
      sprintf(staname, "%s/sta%06d", pathname, ncer);

      // try to open cer* file

      if ( verbose >= VERBOSE_MAXIMAL )
        log( SIGNATURE, "Openning %s\n", cername );

      cerfile.open( cername );

      if ( cerfile.bad() ) 
        error( SIGNATURE, 
               "Cannot open input file: %s\n", cername );

      // @endcode
      // @text
      // Each shower has associated three files: {\em Cherenkov-photons file\/}
      // ({\tt cer*}), {\em Particles file} ({\tt dat*}) and
      // {\em Statistics file} ({\tt sta*}). First we read the data from
      // this last file, which is explained below, and then we read the 
      // {\em Event-Header} block from the {\tt cer*} file.
      // @endtext
      // @code

      // try to open the statistics file

      if ( verbose >= VERBOSE_MAXIMAL )
        log( SIGNATURE, "Openning %s\n", staname );

      stat.openfile ( staname );

      // reading data from this sta* file

      if ( verbose >= VERBOSE_MAXIMAL )
        log( SIGNATURE, "Reading %s\n", staname );

      stat.read();

      // read the Event Header from the cer* file

      evth.read( cerfile );

      // Distinquish between gamma/proton/helium/+heavier nuclei
      if (evth.PrimaryID == gammas) 
            hprimary=0;
      else if (evth.PrimaryID == protons) 
            hprimary=1; 
      else if (evth.PrimaryID == helium) 
            hprimary=2; 
      else
            hprimary=3; 

      // set bin for theta/energy
      
      kbin = ((int)floor(evth.get_theta()/10.)*6 + 
			  (int)floor(log10(evth.get_energy())));

      // get core distance and position

      coreD = evth.get_core(&coreX, &coreY);
          
      t1 = 100000;
      t2 = 0.0;

      ++ntotfiles;

      // initialize time-histogram

      for (i=0; i<60; ++i)
		fstime[kbin][i] = 0.0;

      // read each and every single photon (particle) from cer* file

      rec2[0]  = evth.Etotal;   
      rec2[1]  = coreD;
      rec2[2]  = evth.Theta; 
      rec2[10]  = evth.EvtNumber; 

      if (firstevent == -100) 
        {firstevent = evth.EvtNumber;
         // save random number seed of first event per run
         for (i=0;i<3;++i)
            for (j=0;j<3;++j)
              rnd[i][j]=evth.RndData[i][j];}

      lastevent = evth.EvtNumber; 

      nCphotons = 0;
      while ( ! cerfile.eof() ) {

        // read photon data from the file

        photon.read( cerfile );

        wl = photon.get_wl();

        HF2(300+hprimary,(photon.x-coreX)/100.,(photon.y-coreY)/100.,1.);
        HF1(310+hprimary,sqrt(pow(((photon.x-coreX)/100.),2.)+pow((photon.y-coreY)/100.,2.)),1.);
        HF1(320+hprimary,wl,1.);
        HF1(330+hprimary,photon.h/100000.,1.);

        rec2[3]  = wl; 
        rec2[4]  = photon.x; 
        rec2[5]  = photon.y; 
        rec2[6]  = photon.u; 
        rec2[7]  = photon.v; 
        rec2[8]  = photon.t; 
        rec2[9]  = photon.h; 
//        HFN(hid2,rec2);
	

        if ( wl < 1.0 )
          break;
        
		// fill times

		t = photon.get_t() - stat.get_tfirst();
	
		if ( (t>0.0) && (t<60.0)) {
		  i = (int)t/1.;
		  ++fstime[kbin][i];
		}	  
	
		// hdmytime->Fill(t);
	
        // increase number of Cphotons written

        ++nCphotons;
	
      } // while still there are photons left
      
      // calculate errors for time distribution
      
      for (i=0;i<60;++i)
		fs2time[kbin][i] = sqrt(fstime[kbin][i]);

      // save intermediate vectors in final vector

      for (i=0;i<60;++i) {
		ftstime[kbin][i] += fstime[kbin][i];
		fts2time[kbin][i] += fs2time[kbin][i];
      }
    
      for (i=0;i<NTHSTEP;++i) {
		ftselec[kbin][i] += (stat.plong[2][i] + stat.plong[1][i]);
		fts2elec[kbin][i] += sqrt(stat.plong[2][i] + stat.plong[1][i]);
		ftslight[kbin][i] += stat.plong[8][i];
		fts2light[kbin][i] += sqrt(stat.plong[8][i]);
      }    

      // estimate position of Xmax

      vmax = -1.0;
      for (i=0; i<NTHSTEP; ++i) {
		v = (stat.plong[2][i] + stat.plong[1][i]);  // e- & e+
		if (v > vmax) {
		  vmax = v;
		  j = i;
		}
      }

      // to calculate a rough estimate of Xmax, we make a weigthed mean
      // of the values at maximum, and the previous and next bins
      xmax=-100.;
      if (j > 0) {
         x0 = j*10. + 5.;
         xa = (j-1)*10. + 5.;
         xb = (j+1)*10. + 5.;
         vmaxa = (stat.plong[2][j-1] + stat.plong[1][j-1]); // e- & e+
         vmaxb = (stat.plong[2][j+1] + stat.plong[1][j+1]); // e- & e+
      
         xmax = (xa*vmaxa + x0*vmax + xb*vmaxb) / (vmaxa + vmax + vmaxb);}
      //     else there was no e+e- shower in this event, skip ...

      // fill N-Tuple record

      if (evth.Etotal>0) 
         HF1(100+hprimary,log10(evth.Etotal),1.);
      HF1(110+hprimary,coreD/100.,1.);
      HF1(120+hprimary,evth.Theta/3.14*180.,1.);
      HF1(130+hprimary,evth.Phi/3.14*180.,1.);
      HF2(140+hprimary,coreX/100.,coreY/100.,1.);
      HF1(150+hprimary,evth.zFirstInt/100000.,1.);
      if (stat.get_tlast()-stat.get_tfirst()>0)
         HF1(160+hprimary,stat.get_tlast()-stat.get_tfirst(),1.);
      else 
         HF1(160+hprimary,-5.,1.);

      HF1(170+hprimary,xmax,1.);
      if (evth.Etotal<300){
      if (nCphotons>0) 
         HF1(180+hprimary,log10(nCphotons),1.);
      else
         HF1(180+hprimary,-1.,1.);
      }
      else{
      if (nCphotons>0) 
         HF1(190+hprimary,log10(nCphotons),1.);
      else
         HF1(190+hprimary,-1.,1.);
      }
      rec[0]  = evth.PrimaryID;
      rec[1]  = evth.Etotal;   
      rec[2]  = coreD;
      rec[3]  = evth.Theta; 
      rec[4]  = evth.Phi; 
      rec[5]  = thick(evth.zFirstInt, evth.Theta);   
      rec[6]  = evth.zFirstInt;
      rec[7]  = evth.FirstTarget;
      rec[8]  = evth.p[0];     
      rec[9]  = evth.p[1];     
      rec[10] = evth.p[2];     
      rec[11] = stat.get_tfirst();
      rec[12] = stat.get_tlast();
      rec[13] = xmax;
      rec[14] = coreX;
      rec[15] = coreY;
      rec[16] = nCphotons;

      rec[17] = evth.RunNumber;
      rec[18] = evth.EvtNumber;
      rec[19] = evth.DateRun;
      rec[20] = evth.VersionPGM;
      rec[21] = evth.HeightLev[0];
      rec[22] = evth.SlopeSpec;
      rec[23] = evth.ELowLim; 
      rec[24] = evth.EUppLim; 


      rec[25] = evth.Ecutoffh;  
      rec[26] = evth.Ecutoffm;  
      rec[27] = evth.Ecutoffe;  
      rec[28] = evth.Ecutoffg;  
      
      rec[29] = evth.EGS4yn;  
      rec[30] = evth.NKGyn;
      rec[31] = evth.GHEISHAyn;
      rec[32] = evth.VENUSyn;
      rec[33] = evth.CERENKOVyn;
      rec[34] = evth.NEUTRINOyn;
      rec[35] = evth.HORIZONTyn;
      rec[36] = evth.COMPUTER;
      rec[37] = evth.ThetaMin;
      rec[38] = evth.ThetaMax;
      rec[39] = evth.PhiMin;
      rec[40] = evth.PhiMax;
  
      rec[41] = evth.StepLength;
      rec[42] = evth.CWaveLower;
      rec[43] = evth.CWaveUpper;

      // write record to the file

//      HFN(hid,rec);     
      
      // ntuple->Fill( rec );          
        
      // close files
          
      cerfile.close();
      stat.closefile();
          
      // show how many photons were written
          
      if ( verbose >= VERBOSE_MAXIMAL )
        log( SIGNATURE, "%d C-photons written.\n", nCphotons );

      
    } // while there are still directory entries left

      // check parameters for each RUN...
      // Monitor run statistics:      

      ++rndiag;
      sprintf(ctmp, "%s\n", "-------------------------------------------------------------");
      for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];

      ++rndiag;
      sprintf(ctmp,"%s%f%s%f%s%f\n", "Run: ",evth.RunNumber,
              " Date: ",evth.DateRun, " CorsikaV: ",evth.VersionPGM);
      for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];

      ++rndiag;
      sprintf(ctmp, "%s%s\n", "    Primary Particle: ",
                  primary_name[hprimary]);
      for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];

      ++rndiag;
      sprintf(ctmp, "%s%d%s%d\n", "    First/Last Event Number monitored: ",
                  firstevent," / ",lastevent);
      for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];

      ++rndiag;
      sprintf(ctmp, "%s %f %f %f\n", "    Random number seeds: ",
              rnd[0][0],rnd[0][1],rnd[0][2]);
      for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];

      ++rndiag;
      sprintf(ctmp, "%s %f %f %f\n", "                       : ",
              rnd[1][0],rnd[1][1],rnd[1][2]);
      for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];

      ++rndiag;
      sprintf(ctmp, "%s %f %f %f\n", "                       : ",
              rnd[2][0],rnd[2][1],rnd[2][2]);
      for (i=0;i<256;++i) rdiag[i][rndiag]=ctmp[i];
                        
      // check ranges of i/p params
      //
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                                evth.HeightLev[0],mon_HeightLev, 
                                ": Observation level is "); 
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.SlopeSpec,mon_SlopeSpec[hprimary], 
                    ": Energy-spectrum slope is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.ELowLim,mon_ElowLim[hprimary], 
                    ": Energy-low lim is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.EUppLim,mon_EUppLim[hprimary], 
                    ": Energy-up lim is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.Ecutoffh,mon_Ecutoffh[hprimary], 
                    ": Energy-cutoff hadrons is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.Ecutoffm,mon_Ecutoffm[hprimary], 
                    ": Energy-cutoff muons is ");

      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.Ecutoffe,mon_Ecutoffe[hprimary], 
                    ": Energy-cutoff electrons is ");

      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.Ecutoffg,mon_Ecutoffg[hprimary], 
                    ": Energy-cutoff gammas is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.EGS4yn,mon_EGS4yn[hprimary], 
                    ": EGS4yn flag is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.NKGyn,mon_NKGyn[hprimary], 
                    ": NKGyn flag is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.GHEISHAyn,mon_GHEISHAyn[hprimary], 
                    ": GHEISHAyn flag is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.VENUSyn,mon_VENUSyn[hprimary], 
                    ": VENUSyn flag is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.CERENKOVyn,mon_CERENKOVyn[hprimary], 
                    ": CERENKOVyn flag is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.ThetaMin,mon_ThetaMin[hprimary], 
                    ": ThetaMin  is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.ThetaMax,mon_ThetaMax[hprimary], 
                    ": ThetaMax  is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.PhiMin,mon_PhiMin[hprimary], 
                    ": PhiMin  is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.PhiMax,mon_PhiMax[hprimary], 
                    ": PhiMax  is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.CWaveLower,mon_CWaveLower[hprimary], 
                    ": CWaveLower is ");
      monitor_event(diag,&ndiag,primary_name[hprimary], 
                    evth.CWaveUpper,mon_CWaveUpper[hprimary], 
                    ": CWaveUpper is ");      
    
  } // for each data directory

  // store contents of temporary 1D histograms to profiles

  /*
	Int_t bin;
	Stat_t value;

	for (kbin=0; kbin<18; ++kbin) {

    for (i=0; i<60; ++i) {
	ftstime[kbin][i] /= ntotfiles;
	fts2time[kbin][i] = sqrt(fts2time[kbin][i]/ntotfiles - 
	SQR(ftstime[kbin][i]));
	bin=i+1, value=ftstime[kbin][i];
	cout << bin << ' ' << value << endl << flush;
	htime->SetBinContent( bin, value );
	htime->SetBinError( i+1, fts2time[kbin][i] );
    }
    
    for (i=0; i<NTHSTEP; ++i) {
	ftselec[kbin][i] /= ntotfiles;
	fts2elec[kbin][i] = sqrt(fts2elec[kbin][i]/ntotfiles - 
	SQR(ftselec[kbin][i]));
	helec->SetBinContent( i+1, ftselec[kbin][i] );
	helec->SetBinError( i+1, fts2elec[kbin][i] );
    }

    for (i=0; i<NTHSTEP; ++i) {
	ftslight[kbin][i] /= ntotfiles;
	fts2light[kbin][i] = sqrt(fts2light[kbin][i]/ntotfiles - 
	SQR(ftslight[kbin][i]));
	hlight->SetBinContent( i+1, ftslight[kbin][i] );
	hlight->SetBinError( i+1, fts2light[kbin][i] );
    }

    TH1F h1;
    TH1F h2;
    TH1F h3;

    htime->Copy( h1 );
    helec->Copy( h2 );
    hlight->Copy( h3 );

    hfile.Write();
	
	}
  */
  
  // close output file
  
  if ( verbose >= VERBOSE_MINIMAL ) 
	log( SIGNATURE, "Closing N-Tuple file %s\n", outname );
  
  // hfile.Write();
  // hfile.Close();  
  
  HROUT(0,icycle," ");

  fprintf(stdout, "%s\n", "-------------------------------------------------------------");
  fprintf(stdout, "%s\n", "---------------------- RUN STATISTICS -----------------------");

  for (i=0;i<=rndiag;++i) 
      {char tmp[256]; int j;
      {for (j=0;j<256;++j) tmp[j]=rdiag[j][i];}
       fprintf(stdout,"%s",tmp);}

  fprintf(stdout, "%s\n", "-------------------------------------------------------------");
  fprintf(stdout, "%s\n", "----------------------- DIAGNOSTICS -------------------------");
  fprintf(stdout, "%s\n", "-------------------------------------------------------------");

  if (ndiag<0) 
     {fprintf(stdout, "%s\n", " None, everything's splendidly perfect ...");
      fprintf(stdout, "%s\n", " Nada, todo fue bien, hacemos fiesta ...");
      fprintf(stdout, "%s\n", " Keine Fehler gefunden ...");}
  else 
    {for (i=0;i<=ndiag;++i) 
         {char tmp[256]; int j;
         {for (j=0;j<256;++j) tmp[j]=diag[j][i];}
          fprintf(stdout,"%s",tmp);}}


  HREND("CHECKMC");
  
  // program finished
  
  if ( verbose >= VERBOSE_MINIMAL ) 
  	log( SIGNATURE, "Done.\n");
  
}
// @endcode

// @T \newpage

// @subsection Functions definition

//------------------------------------------------------------
// @name present
//
// @desc Make some presentation
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//------------------------------------------------------------
// present
//
// make some presentation
//------------------------------------------------------------

void 
present(void)
{
cout << "##################################################\n"
<<  SIGNATURE << '\n' << '\n'
<< "Check of the MC data generated by MMCS\n"
<< "J C Gonzalez, P Morawitz, D Petry, 27 Oct 1998\n"
<< "##################################################\n\n";
}
// @endcode

//------------------------------------------------------------
// @name log                             
//                                   
// @desc function to send log information
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//------------------------------------------------------------
// log
//
// function to send log information
//------------------------------------------------------------

void
log(const char *funct, char *fmt, ...)
{
  va_list args;
  
  //  Display the name of the function that called error
  printf("[%s]: ", funct);
  
  // Display the remainder of the message
  va_start(args, fmt);
  vprintf(fmt, args);
  va_end(args);
}
// @endcode

//------------------------------------------------------------
// @name error                                                    
//                                                           
// @desc function to send an error message, and abort the program
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//------------------------------------------------------------
// error
//
// function to send an error message, and abort the program
//------------------------------------------------------------

void
error(const char *funct, char *fmt, ...)
{
  va_list args;

  //  Display the name of the function that called error
  fprintf(stderr, "ERROR in %s: ", funct);

  // Display the remainder of the message
  va_start(args, fmt);
  vfprintf(stderr, fmt, args);
  va_end(args);

  fprintf(stderr, "\nTerminating program.\n");
  exit(999);
}
// @endcode


//------------------------------------------------------------
// @name readparam                                        
//                                                    
// @desc function to read the parameters of the check
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//------------------------------------------------------------
// readparam
//
// function to read the parameters of the check
//------------------------------------------------------------

void 
readparam(void)
{
  char sign[] = GLUE_postp( PROGRAM, VERSION ); // initialize sign
  char line[LINE_MAX_LENGTH];    // line to get from the stdin
  char token[ITEM_MAX_LENGTH];   // a single token
  int i, j;                      // dummy counters

  // get signature
  cout << "\nExpecting input from stdin (see example input file) ...\n";
  cin.getline(sign, strlen(SIGNATURE) + 1);
  //  cin.getline(line, LINE_MAX_LENGTH);
  if (strcmp(sign, SIGNATURE) != 0) {
    cerr << "ERROR: Signature of parameters file is not correct\n";
    cerr << '"' << sign << '"' << '\n';
    cerr << "should be: " << SIGNATURE << '\n';
    exit(1);
  }

  // loop till the "end" directive is reached
  int is_end = FALSE;
  while (! is_end) {          

    // get line from stdin
    cin.getline(line, LINE_MAX_LENGTH);

    // look for each item at the beginning of the line
    for (i=0; i<=end_file; i++) 
      if (strstr(line, PAR_ITEM_NAMES[i]) == line)
        break;
        
    // if it is not a valid line, just ignore it
    if (i == end_file+1) {
      cerr << "    Skipping unknown token in [" << line << "]\n";
      continue;
    }

    // case block for each directive
    switch ( i ) {

    case data_paths:          // beginning of the list of valid paths
          
      //get number of paths to read
      sscanf(line, "%s %d", token, &Num_of_paths);

      if ( verbose == TRUE )
        cout << "    " << Num_of_paths << " path(s) will be used.\n";

      // allocate memory for paths list
      Paths_list = (char**)malloc(Num_of_paths * sizeof(char*));
      for (i=0; i<Num_of_paths; i++) 
        Paths_list[i] = (char*)malloc(PATH_MAX_LENGTH * sizeof(char));

      // read each single path
      for (i=0; i<Num_of_paths; i++) {

        cin.getline(Paths_list[i], PATH_MAX_LENGTH);

        // remove empty spaces at the end
        for (j=0; j<strlen(Paths_list[i]); j++) {
          if (Paths_list[i][j] == ' ') {
            Paths_list[i][j] = '\0';
            break;
          }                     
        }

        // show the path
        if ( verbose == TRUE )
          cout << "    " << '[' << Paths_list[i] << ']' << '\n';
      }

      break;

    case output_file:         // name of the output file
          
      // get the name of the output_file from the line
      sscanf(line, "%s %s", token, Output_filename);

      break;     

    case verbose_level:      // verbose level 0-4
          
      // get the verbose level
      sscanf(line, "%s %d", token, &verbose);
          
      break;

    case max_events:          // maximum number of events
          
      // get the number of events
      sscanf(line, "%s %d", token, &max_Events);
          
      break;

    case end_file:            // end of the parameters file

      // show a short message
      is_end = TRUE;

      break;

    } // switch ( i ) 

  } // while (! is_end)

  // after the loop is finished, return to the main function
  return;
}
// @endcode


//------------------------------------------------------------
// @name get_num_of_paths                 
//                                                
// @desc get number of defined paths
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//--------------------------------------------------
// get_num_of_paths
//
// get number of defined paths
//--------------------------------------------------
int
get_num_of_paths(void)
{
  return (Num_of_paths);
}
// @endcode


//------------------------------------------------------------
// @name get_path_name
//                                                
// @desc get path name number "i" from the list of paths
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//--------------------------------------------------
// get_path_name
//
// get path name number "i" from the list of paths
//--------------------------------------------------
char *
get_path_name(int i)
{
  return (Paths_list[i]);
}
// @endcode


//------------------------------------------------------------
// @name get_output_filename
//                                                
// @desc get name of the output file
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//--------------------------------------------------
// get_output_filename
//
// get name of the output file
//--------------------------------------------------
char *
get_output_filename(void)
{
  return (Output_filename);
}
// @endcode




//------------------------------------------------------------
// @name get_verbose
//                                                
// @desc get the "verbose" status
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//--------------------------------------------------
// get_verbose
//
// get the "verbose" status
//--------------------------------------------------
int
get_verbose(void)
{
  return (verbose);
}
// @endcode


//------------------------------------------------------------
// @name get_max_event              
//                                          
// @desc get max events number
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//--------------------------------------------------
// get_max_event
//
// get max events number
//--------------------------------------------------
int 
get_max_events(void)
{
  return( max_Events );
}
// @endcode


//------------------------------------------------------------
// @name thick       
//                                          
// @desc converts height [cm] -> thickness [g cm-2]
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//--------------------------------------------------
// thick
//
// converts height [cm] -> thickness [g cm-2]
//--------------------------------------------------

float 
thick( float height, float theta )
{
  float h, thk;

  // Ibarra's modification
  h = sqrt( SQR(REarth) + 
            SQR(height/cos(theta)) + 
            (2.0*REarth*height)) - REarth;

  if     ( h < 4.e5 ) {
    thk = AATM[0] + BATM[0] * exp ( -h * DATM[0] );
  } else if ( h < 1.e6 ) {
    thk = AATM[1] + BATM[1] * exp ( -h * DATM[1] );
  } else if ( h < 4.e6 ) {
    thk = AATM[2] + BATM[2] * exp ( -h * DATM[2] );
  } else if ( h < 1.e7 ) {
    thk = AATM[3] + BATM[3] * exp ( -h * DATM[3] );
  } else {
    thk = AATM[4] - h * CATM[4];
  }

  return ( thk );

}
// @endcode

//------------------------------------------------------------
// @name height   
//                                          
// @desc converts thickness [g cm-2] -> height [cm]
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//--------------------------------------------------
// height
//
// converts thickness [g cm-2] -> height [cm]
//--------------------------------------------------

float 
height( float thk, float theta )
{
  float h, height;

  if     ( thk > 631.1e0 ) {
    h = CATM[0] * log ( BATM[0] / (thk - AATM[0]) );
  } else if ( thk > 271.7e0 ) {
    h = CATM[1] * log ( BATM[1] / (thk - AATM[1]) );
  } else if ( thk > 3.0395e0 ) {
    h = CATM[2] * log ( BATM[2] / (thk - AATM[2]) );
  } else if ( thk > 0.00128292e0 ) {
    h = CATM[3] * log ( BATM[3] / (thk - AATM[3]) );
  } else {
    h = (AATM[4] - thk) * DATM[4];
  }
  
  // Ibarra's modification
  height = SQR(cos(theta)) * 
    (-REarth + sqrt(SQR(REarth) + 
                    (SQR(h + (2.0*REarth*h))/SQR(cos(theta)))));

  return ( height );

}
// @endcode
 
//------------------------------------------------------------
// @name monitor_event   
//                                          
// @desc compare parameter with reference and report if not equal
//
// @date Sat Jun 27 05:58:56 MET DST 1998
// @function @code 
//--------------------------------------------------
// monitor_event
//
// compare parameter with reference
//--------------------------------------------------
void monitor_event(char diag[DIAGX][DIAGY], int *ndiag, const char* primary, 
                          const float a, const float b, const char* message) 
{char diagt[DIAGX]; int i; 
 if (a!=b){
   if (*ndiag<(DIAGY-1)) 
     {++*ndiag;
      sprintf(diagt,"%s%s%s%s%f%s%f\n","Primary ",primary," ",
	      message,a," but should be ",b);} 
   else{
     sprintf(diagt,"%s",
	"Two many errors encountered, won't print any further diagnostics\n");
     for (i=0;i<DIAGX;++i){ 
       diag[i][*ndiag]=diagt[i];
     }
   } // endif
 } // endif
 return;
}