/*
 *MM        MM UU        UU PPPPPPPPPPP   AAAAAAAAAA   GGGGGGGGGG  EEEEEEEEEEEE
 *MMM      MMM UU        UU PPPPPPPPPPPP AAAAAAAAAAAA GGGGGGGGGGGG EEEEEEEEEEEE
 *MMMM    MMMM UU        UU PP        PP AA        AA GG        GG EE
 *MM MM  MM MM UU        UU PP        PP AA        AA GG           EE
 *MM  MMMM  MM UU        UU PP        PP AA        AA GG           EE
 *MM   MM   MM UU        UU PPPPPPPPPPPP AAAAAAAAAAAA GG           EEEEEEEE
 *MM        MM UU        UU PPPPPPPPPPP  AAAAAAAAAAAA GG     GGGGG EEEEEEEE
 *MM        MM UU        UU PP           AA        AA GG     GGGGG EE
 *MM        MM UU        UU PP           AA        AA GG        GG EE
 *MM        MM UU        UU PP           AA        AA GG        GG EE
 *MM        MM UUUUUUUUUUUU PP           AA        AA GGGGGGGGGGGG EEEEEEEEEEEE
 *MM        MM  UUUUUUUUUU  PP           AA        AA  GGGGGGGGGG  EEEEEEEEEEEE
 *
 *@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
 * 
 * Author:  G.Carminati
 * Version: v3r5 - 16th April 2010
 * 
 */


#include<iostream>
#include<fstream>
#include<cmath>
#include<iomanip>
#include<string>
#include<sstream>

#include "Parameters.hh"
#include "Date.hh"
#include "Muons.hh"
#include "Decode.hh"
#include "OutputWriter.hh"

#include "TMath.h"
#include "TRandom3.h"

#ifdef _ANTARES_ENABLED__
#include "io_gcc.hh"          /* by J.Brunner */
#endif
#ifdef _ANTCC_ENABLED__
#include "ReadDetector.hh"
#endif

void PrintHeader(void);

static const double c_light_ns = TMath::C()* 1.e-9;  	/* in m/ns */

int main(int argc, char** argv)
{

  /* print code header on screen */
  PrintHeader();
  
  /* information for execution */
  lib_det = false;
  histo   = false;
  Decode dec;
  int nerr = dec.CommandLine(argc, argv);
  if(nerr < 0)
    {
      std::cout << dec.ErrorCommandLine();
      exit(3);
    }
  int nerror = 0;
  nerror = dec.DecodeParameters();
  if(nerror > 0)
    {
      std::cout << dec.ErrorDecoder(nerror);
      exit(2);
    }
  else if(nerror < 0)
    {
      std::cout << dec.ErrorOpenFiles(nerror);
      exit(1);
    }

  
  OutputWriter OutWrt;  


  /* opening output files */
  std::ofstream livetime_fout;
  
  livetime_fout.open(LivetimeFileName.c_str(), std::ios::out);

/*  to be moved in the class?
  if (fout.fail() )
    {
      nerror = -98; 
      std::cout << dec.ErrorOpenFiles(nerror);
      exit(1);
    }
  if (livetime_fout.fail() )
    {
      nerror = -97; 
      std::cout << dec.ErrorOpenFiles(nerror);
      exit(1);
    }
*/


  /* seeds definition */
  if(seed == 0)
    seed = time(NULL);
  TRandom3 genrandom(seed);

  /* Load geometry data from the detector file */
  if (lib_det)
    {
#ifdef _ANTCC_ENABLED__
      ReadDetector detector(DetectorFileName);
      DEPTHmax = detector.surface_z * MeterToKm();
      Vec3D cog = detector.compute_cg();

      double z_min = 0.;
      double z_max = 0.;
      Vec2D z_det(z_min, z_max);
      double   r_det = 0.;
      double  OffSet = 0.;
      double   h_can = 0.;
      double   r_can = 0.;

      /* read from detector file the dimension of detector */
      r_det =  detector.compute_radius(cog);
      z_det = detector.compute_Z(cog);

      /* compute from detector file the minimum and the maximum value 
	 of CAN z-coordinate and compute the CAN radius */ 
      OffSet = fabs(NAbsLength*AbsLength);
      z_min = max(detector.ground_z-cog.z, z_det.x-OffSet);
      z_max = min(detector.surface_z-cog.z, z_det.y+OffSet);
      h_can = z_max - z_min;
      r_can = r_det + OffSet;

      Zmin = z_min;
      Zmax = z_max;
      CANr = r_can;
      CANh = h_can;
#endif 
    }

  /* create parameters for the CAN */
  Vec3D can(Zmin, Zmax, CANr);
  if(!lib_det)
    CANh = CANheight();

  /* the CAN have to be enlarged to avoid losing events */
  Vec3D enlargedCAN(can.x, can.y, can.z+EnlargedCANr);

  /* check on vertical depth values */
  if (weDEPTHmin() < 1.5 || weDEPTHmin() > 5.)
    {
      nerror = 6;
      std::cout << dec.ErrorDecoder(nerror);
      exit(2);
    }
  if (weDEPTHmax() < 1.5 || weDEPTHmax() > 5.)
    {
      nerror = 7;
      std::cout << dec.ErrorDecoder(nerror);
      exit(2);
    }


  OutWrt.WriteHeader();  

  /* Implementation class Muons */
  Muons muon (Qmin(), 0., MULTmin, weDEPTHmin(), genrandom);
  muon.GetParameters(DEPTHmax, enlargedCAN.x, enlargedCAN.y, 
		     enlargedCAN.z, Rmin, Rmax, Emin, Emax);

  /* initialization of the counter to compute the generated events 
     on the enlarged can */
  unsigned int counter = 0;

  /* initialization of counters to compute livetime */
  unsigned int thetarange = 86;
  unsigned int ncounter[thetarange];
  memset(ncounter, 0, thetarange*sizeof(int));

  /* compute the area of the enlarged CAN (in m**2) */
  double Sz  = TMath::Pi()*enlargedCAN.z*enlargedCAN.z;      	//<------ upper disk
  double Sxy = 2*enlargedCAN.z*CANh;                		//<------ barrel surface

  /* compute the maximum projected area using the fact that 
     d(Area)/d(theta) = 0 (in m**2) */
  double ALFA = atan( - enlargedCAN.z /( 2*CANh ))*.5 + TMath::Pi()/2. ;

  /* compute the maximum number of events per unit of time */
  double halfdegree = 0.5 * DegToRad();
  double Smax = Sz*cos(ALFA) + Sxy*sin(ALFA);
  Smax = max(Sz, Smax);
  double Nmax = muon.Flux()*Smax*muon.SolidAngle(ALFA-halfdegree,
						 ALFA+halfdegree);

  /* counter for events with energy lower than threshold one */
  unsigned int Nlow = 0;

  for (unsigned int i = 0; i < events; ++i)
    {
      muon.GetParameters(DEPTHmax, enlargedCAN.x, enlargedCAN.y, 
			 enlargedCAN.z, Rmin, Rmax, Emin, Emax);

      unsigned int m = 0;
      double phi     = 0.;
      /* thetaSG means the zenith angle in the System of Generation */
      double thetaSG = 0.;
      double ctSG    = 0.;
      double stSG    = 0.;
      /* thetaDF means the zenith angle in the Direction-to-fly System */
      double thetaDF = 0.;
      double ctDF    = 0.;
      double stDF    = 0.;
      double u       = 0.;
      double wedepthBundle = 0.;
      Vec3D posBundle(0.,0.,0.);
      double evttime = 0.;

      /* generate random the azimuth (phi) angle [in radians] */
      phi = genrandom.Rndm() * 2.*TMath::Pi();
      /* WARNING: the azimuthal angle is the angle of the pointlike
       *          source, so upward-going; instead, the muon is
       *          downward-going. So, it is necessary to convert
       *          the azimuthal angle into
       *          phi_mu = 180 degrees + azimuth
       */
      phi += TMath::Pi();    /* in rad */
      
      /* generate the event time is requested */
      
      if(GenTime)
        evttime = double(TimeStampStart) + (double)(TimeStampStop - TimeStampStart) * genrandom.Rndm();

      bool loop = true;
      while(loop)
	{
	  /* generate random the muon multiplicity */
	  m = (int)(genrandom.Rndm()*(MULTmax + 1 - MULTmin)) + MULTmin;

	  /* generate random the zenith (theta) angle [in radians] */
 	  thetaSG = genrandom.Rndm() * ( Qmax() - Qmin() ) + Qmin();
	  ctSG = cos(thetaSG);
	  stSG = sin(thetaSG);

	  /* implement new parameters in the class */
	  muon.Set(thetaSG, phi, m, wedepthBundle);

	  /* compute the shower axis position 
	     on the can surface (in m) */
	  posBundle = muon.axisBundleOnCAN();

	  /* convert the z-position of the impact point 
	     in the vertical depth with respect to the sea level 
	     (in km.w.e.) */
	  wedepthBundle = muon.Depth(posBundle.z);

	  /* compute the projected area of the enlarged can 
	     for the generated zenithal angle (in m**2) */
 	  double S = Sz*ctSG + Sxy*stSG;

	  /* compute the number of events per unit of time 
	     for generated zenithal angle, multiplicity and vertical depth */
 	  double Nproj = muon.Flux()*S*muon.SolidAngle(thetaSG+halfdegree,
						       thetaSG-halfdegree);

	  /* accept or reject event */
	  u = genrandom.Rndm() * Nmax;
	  if (u <= Nproj)         //<------- accept event
	    loop = false;
	}

      /* compute the generated single events N* 
	 on upper disk of the enlarged CAN at different zenithal angles */
      float fhBundle = wedepthBundle;
      float    fHmin = weDEPTHmin();
      if (fhBundle == fHmin && m == MULTmin)
	{
	  ++counter;
	  for(unsigned int j = 0; j < thetarange; ++j)
	    {
	      if(thetaSG * RadToDeg() > j && thetaSG * RadToDeg() <= j+1)
		++(ncounter[j]);
	    }
	}

      /* WARNING: Accordling with the switching of the azimuthal angle, 
       *          the zenithal angle must be switched into
       *          theta_mu = 180 degrees - zenith
       */
      thetaDF = TMath::Pi() - thetaSG;     /* in rad */
      ctDF = cos(thetaDF);
      stDF = sin(thetaDF);
      muon.Set(thetaSG, phi, m, wedepthBundle);


      /* unit vector (i.e. direction cosines) of Bundle Axis */
      double cf = cos(phi);
      double sf = sin(phi);
      Vec3D uBundle(stDF*cf, stDF*sf, ctDF);

      if (m == 1) /* single muons */
      	{
	  double t = 0.;
	  int MuOnCAN=0;
		      
	  /* compute the energy of single muon */
	  double e = muon.HoR_Energy_singlemu();    /* in TeV */

	  /* check if the generated energy is greater than the 
	     energy of threshold */
	  if (e >= Ethreshold)
	    {
	      OutWrt.WriteEvent(i, m, posBundle, uBundle, e, &posBundle.x, &posBundle.y, &posBundle.z, &e, &t, &MuOnCAN, evttime); 	      
	    }
	  else /* events under the energy of threshold */
	    {
	      ++Nlow;
	      --i;
	      continue;
	    }
	}

      else /* multiple muons */

      	{
	  unsigned int mc = 0;
	  double bundle_energy = 0.;
	  double bundle_energy_at_can = 0.; 
	  double xMuOnLAB[m];
	  double yMuOnLAB[m];
	  double zMuOnLAB[m];
	  double xMuOnCAN[m];
	  double yMuOnCAN[m];
	  double zMuOnCAN[m];
	  double distance[m];
	  double delaytime[m];
	  double E[m];	  
	  int MuOnCAN[m];  
	  
	  memset(xMuOnLAB, 0, m*sizeof(double));
	  memset(yMuOnLAB, 0, m*sizeof(double));
	  memset(zMuOnLAB, 0, m*sizeof(double));
	  memset(xMuOnCAN, 0, m*sizeof(double));
	  memset(yMuOnCAN, 0, m*sizeof(double));
	  memset(zMuOnCAN, 0, m*sizeof(double));
	  memset(distance, 0, m*sizeof(double));
	  memset(delaytime, 0, m*sizeof(double));
	  memset(E, 0, m*sizeof(double));
	  memset(MuOnCAN, 0, m*sizeof(int));

      	  for (unsigned int l = 0; l < m; ++l)
      	    {
	      /* compute the lateral spread of each muon in bundle */
      	      double R[m];
	      memset(R, 0, m*sizeof(double));
	      R[l]= muon.HoR_Radius();

	      /* compute the position of each muon in bundle on the laboratory 
		 frame (in m) */
 	      Vec3D posMuOnLAB = muon.multimuOnLAB(posBundle, R[l]);
	      xMuOnLAB[l] = posMuOnLAB.x;
	      yMuOnLAB[l] = posMuOnLAB.y;
	      zMuOnLAB[l] = posMuOnLAB.z;

 	      /* compute if the i-th muon in bundle impacts point also
		 on the can surface */ 
	      Vec3D posMuOnCAN (0.,0.,0.);
	      int iok = muon.muOnCAN(posMuOnLAB, uBundle, posMuOnCAN);
	      MuOnCAN[l] = iok;
	      //if (iok != 0 && l == 0) // muon does not intercept the can
                //std::cout << "this muon would be forced to enter the can" << std::endl;
		//--l;   V.Kulikovskiy - why would we need to force the first muon in the bundle intersect the can
	      //else 
              if (iok == 0) /* muon intercepts the can */
		{
		  xMuOnCAN[l] = posMuOnCAN.x;
		  yMuOnCAN[l] = posMuOnCAN.y;
		  zMuOnCAN[l] = posMuOnCAN.z;
		
		  /* compute the delay time (in ns) */
		  distance[l]=sqrt( (xMuOnCAN[l] - xMuOnLAB[l])*
				    (xMuOnCAN[l] - xMuOnLAB[l]) +
				    (yMuOnCAN[l] - yMuOnLAB[l])*
				    (yMuOnCAN[l] - yMuOnLAB[l]) +
				    (zMuOnCAN[l] - zMuOnLAB[l])*
				    (zMuOnCAN[l] - zMuOnLAB[l]) );
		  if (zMuOnLAB[l] < zMuOnCAN[l])
		    distance[l] = - distance[l];
                  //V. Kulikovskiy. It is better to have time respect to the bundle time on the can so the further should be remove. Also, note that it is not valid for the events without fist muon reaching the can (and that is why it was forced to reach it!)
/*
		  if (l == 0)
		    delaytime[l] = 0.;
		  else
		    delaytime[l] = (distance[l] - distance[0])/c_light_ns;
*/

		  delaytime[l] = distance[l]/c_light_ns; 	/* in ns */
		      
		  /* compute the energy of multiple muons */
		  // BEGIN approach B
		  E[l] = muon.HoR_Energy_multimu(R[l]);
                  bundle_energy += E[l];

		  //E[l] -= distance[l]*0.2*0.001;  //V. Kulikovskiy - very rough energy correction assuming 2MeV/cm (0.2GeV/m)
		  //if (E[l] < 0) E[l] = 0;
		  E[l] = muon.calcEloss(E[l],distance[l]);
		  bundle_energy_at_can += E[l];
		  ++mc;
		  
		} else {
                //V. Kulikovskiy muon bundle energy should contain total energy including muons not intercecting the can
                E[l] = muon.HoR_Energy_multimu(R[l]);
                bundle_energy += E[l];
                }
	    }


	  if (bundle_energy_at_can >= Ethreshold){
	    OutWrt.WriteEvent(i, m, posBundle, uBundle, bundle_energy, xMuOnCAN, yMuOnCAN, zMuOnCAN, E, delaytime, MuOnCAN, evttime);  
	  } 
	  else /* events under the energy of threshold */
	    {
	      ++Nlow;
	      --i;
	      continue;
	    }
	}   
    }

  /* compute livetime */
  livetime_fout << std::endl;
  livetime_fout <<"Generated events: " << events << std::endl;
  livetime_fout <<"Generated events with multiplicity " << MULTmin 
		<< " on the upper disk of the CAN: " << counter << std::endl;
  livetime_fout << "\nEvents with energy higher than energy of threshold: "
		<< events <<std::endl;
  livetime_fout << "Events with energy lower than energy of threshold: " 
		<< Nlow << std::endl;
  livetime_fout << std::endl;
  livetime_fout << "theta (degrees)    N*   Delta N*    N*_MC     T (sec) "
		<< "  Delta T\n";
  livetime_fout << "******************************************************"
		<< "**********\n";

  /* range of 0 to 9 degrees */
  unsigned int nsmallangle = 0;
  for(unsigned int k = 0; k < 10; ++k)
    nsmallangle += ncounter[k];
  unsigned int ensmallangle = (int)(sqrt((double)nsmallangle)+0.5);
  muon.Set(5*DegToRad(), 0., MULTmin, weDEPTHmin());
  muon.GetParameters(DEPTHmax, enlargedCAN.x, enlargedCAN.y, 
		     enlargedCAN.z, Rmin, Rmax, Emin, Emax);
  double N_MC = muon.GAMMA(0.*DegToRad(),10.*DegToRad());
  double livetime10  = (double)nsmallangle/N_MC;
  double elivetime10 = (double)ensmallangle/N_MC;
  double weight10 = 0.;
  if(livetime10 != 0)
    weight10 = 1./(elivetime10*elivetime10);
  livetime_fout.setf(std::ios::fixed|std::ios::right);
  livetime_fout	<< std::setw(11) << "[0,10]" << std::setw(10) << nsmallangle 
		<< std::setw(9) << ensmallangle 
		<< std::setprecision(1) << std::setw(10) << N_MC;
  livetime_fout.unsetf(std::ios::fixed);
  livetime_fout.setf(std::ios::floatfield); 
  livetime_fout <<std::setprecision(3)<<std::setw(11)<<livetime10 
		<<std::setprecision(1)<<std::setw(10)<<elivetime10<<std::endl;

  /* range of 10 to 69 degrees */
  unsigned int  twocounter = 0;
  unsigned int etwocounter = 0;
  double livetime[thetarange];
  double elivetime[thetarange];
  double weight[thetarange];
  memset(livetime, 0, thetarange*sizeof(double));
  memset(elivetime, 0, thetarange*sizeof(double));
  memset(weight, 0, thetarange*sizeof(double));

  for(unsigned int jj = 10; jj < 70; ++jj)
    {
      double FASA = 0.;
      if(jj % 2 != 0)
	{
	  twocounter = ncounter[jj-1] + ncounter[jj];
	  etwocounter = (int)(sqrt((double)twocounter)+0.5);
	  muon.Set( (double)(jj*DegToRad()), 0., MULTmin, weDEPTHmin());
	  FASA = muon.GAMMA( (double)(jj-1)*DegToRad(),
			     (double)(jj+1)*DegToRad());
	  livetime[jj]  = (double)twocounter/FASA;
	  elivetime[jj] = (double)etwocounter/FASA;
	  weight[jj] = 1./(elivetime[jj]*elivetime[jj]);
	  livetime_fout.unsetf(std::ios::floatfield);
	  livetime_fout.setf(std::ios::fixed|std::ios::right);
	  livetime_fout << std::setw(5) << "]" <<jj-1 << ","<< jj+1 <<"]"
			<< std::setw(10) << twocounter 
			<< std::setw(9) << etwocounter 
			<< std::setprecision(1) << std::setw(10) << FASA;
	  livetime_fout.unsetf(std::ios::fixed);
	  livetime_fout.setf(std::ios::floatfield);
	  livetime_fout <<std::setprecision(3)<<std::setw(11)<<livetime[jj] 
			<<std::setprecision(1)<<std::setw(10)<<elivetime[jj] 
			<<std::endl;
	}
    }

  /* range of 70 to 75 degrees */
  unsigned int nangle70 = 0;
  for(unsigned int kk = 70; kk < 76; ++kk)
    nangle70 += ncounter[kk];
  unsigned int enangle70 = (int)(sqrt((double)nangle70)+0.5);
  muon.Set(72.5*DegToRad(), 0., MULTmin, weDEPTHmin());
  double N_MC70 = muon.GAMMA(70.*DegToRad(),75.*DegToRad());
  double livetime70  = (double)nangle70/N_MC70;
  double elivetime70 = (double)enangle70/N_MC70;
  double weight70 = 0.;
  if(livetime70 != 0)
    weight70 = 1./(elivetime70*elivetime70);
  livetime_fout.unsetf(std::ios::floatfield);
  livetime_fout.setf(std::ios::fixed|std::ios::right);
  livetime_fout << std::setw(11) << "]70,76]" << std::setw(10) << nangle70 
		<< std::setw(9) << enangle70 
		<< std::setprecision(1) << std::setw(10) << N_MC70;
  livetime_fout.unsetf(std::ios::fixed);
  livetime_fout.setf(std::ios::floatfield);
  livetime_fout	<<std::setprecision(3)<<std::setw(11)<<livetime70
		<<std::setprecision(1)<<std::setw(10)<<elivetime70<<std::endl;

  /* range of 76 to 85 degrees */
  unsigned int nangle76 = 0;
  for(unsigned int k76 = 76; k76 < 86; ++k76)
    nangle76 += ncounter[k76];
  unsigned int enangle76 = (int)(sqrt((double)nangle76)+0.5);
  muon.Set(80.*DegToRad(), 0., MULTmin, weDEPTHmin());
  double N_MC80 = muon.GAMMA(76.*DegToRad(),85.*DegToRad());
  double livetime80  = (double)nangle76/N_MC80;
  double elivetime80 = (double)enangle76/N_MC80;
  double weight80 = 0.;
  if(livetime80 != 0)
    weight80 = 1./(elivetime80*elivetime80);
  livetime_fout.unsetf(std::ios::floatfield);
  livetime_fout.setf(std::ios::fixed|std::ios::right);
  livetime_fout << std::setw(11) << "]76,85]" << std::setw(10) << nangle76 
		<< std::setw(9) << enangle76 
		<< std::setprecision(1) << std::setw(10) << N_MC80; 
  livetime_fout.unsetf(std::ios::fixed);
  livetime_fout.setf(std::ios::floatfield);
  livetime_fout	<<std::setprecision(3)<<std::setw(11)<<livetime80
		<<std::setprecision(1)<<std::setw(10)<<elivetime80<<std::endl; 

  double sumweight = 0.;
  double num = 0.;
  double sumnum = 0.;
  double avlivetime = 0.;
  double eavlivetime = 0.;
  for(unsigned int hh = 10; hh < 70; ++hh)
    {
      if(ncounter[hh] != 0)
	{
	  num = livetime[hh] * weight[hh];
	  sumnum += num;
	  sumweight += weight[hh];
	}
    }
  sumnum += livetime10*weight10 + livetime70*weight70 + livetime80*weight80;
  sumweight += weight10 + weight70 + weight80;

  avlivetime = sumnum/sumweight;
  eavlivetime = 1./(sqrt(sumweight));

  livetime_fout.unsetf(std::ios::scientific);
  livetime_fout.setf(std::ios::fixed|std::ios::floatfield);
  livetime_fout << std::endl;
  livetime_fout << std::setprecision(3) << "Average livetime = (" <<avlivetime 
		<< " pm " << std::setprecision(1) << eavlivetime << ") sec\n";
  livetime_fout << std::setprecision(3) << "                 = (" 
		<< avlivetime*sTOdays() << " pm " << std::setprecision(1)
		<< eavlivetime*sTOdays() << ") days\n";
  livetime_fout << std::endl;


  OutWrt.WriteLiveTime(avlivetime,eavlivetime);

  livetime_fout.close();
}

///////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void PrintHeader(void)
{

std::cout << std::endl;
std::cout << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;
std::cout << "                                                                                      " << std::endl;
std::cout << "                                                                                      " << std::endl;
std::cout << "     MM        MM UU        UU PPPPPPPPPPP   AAAAAAAAAA   GGGGGGGGGG  EEEEEEEEEEEE    " << std::endl;
std::cout << "     MMM      MMM UU        UU PPPPPPPPPPPP AAAAAAAAAAAA GGGGGGGGGGGG EEEEEEEEEEEE    " << std::endl;
std::cout << "     MMMM    MMMM UU        UU PP        PP AA        AA GG        GG EE              " << std::endl;
std::cout << "     MM MM  MM MM UU        UU PP        PP AA        AA GG           EE              " << std::endl;
std::cout << "     MM  MMMM  MM UU        UU PP        PP AA        AA GG           EE              " << std::endl;
std::cout << "     MM   MM   MM UU        UU PPPPPPPPPPPP AAAAAAAAAAAA GG           EEEEEEEE        " << std::endl;
std::cout << "     MM        MM UU        UU PPPPPPPPPPP  AAAAAAAAAAAA GG     GGGGG EEEEEEEE        " << std::endl;
std::cout << "     MM        MM UU        UU PP           AA        AA GG     GGGGG EE              " << std::endl;
std::cout << "     MM        MM UU        UU PP           AA        AA GG        GG EE              " << std::endl;
std::cout << "     MM        MM UU        UU PP           AA        AA GG        GG EE              " << std::endl;
std::cout << "     MM        MM UUUUUUUUUUUU PP           AA        AA GGGGGGGGGGGG EEEEEEEEEEEE    " << std::endl;
std::cout << "     MM        MM  UUUUUUUUUU  PP           AA        AA  GGGGGGGGGG  EEEEEEEEEEEE    " << std::endl;
std::cout << "                                                                                      " << std::endl;
std::cout << "                                                                                      " << std::endl;
std::cout << "                                   "<<versn<<"                                        " << std::endl;
std::cout << "                                                                                      " << std::endl;
std::cout << "             the fast atmospheric MUon GEnerator for neutrino telescopes              " << std::endl;
std::cout << "                           based on PArametric formulas                               " << std::endl;
std::cout << "                                                                                      " << std::endl;
std::cout << "                                                                                      " << std::endl;
std::cout << " @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@" << std::endl;


  return;
}

/*! The MUPAGE code was developed for the production of atmospheric
multimuons on the underwater detector CAN surface, without a
full Monte Carlo simulation but using parametric formulas.
MUPAGE is for ANTARES users.
*/