#ifndef CDASUNIV_H
#define CDASUNIV_H

#include "UnivParam.h"

// Usage: UnivCDAS().GetSignal(...) returns signal in VEM, see parameters below

class CDASUniv {
  public:
    int atmos;         // atmos: month of the year, -1: average
    double rhoGround;     // rhoGround: density at ground
    double GetSignal(
      double logE, 
      double Xmax,
      double Nmu, 
      double theta, 
      double z,
      double r,
      double psi, 
      int component=-1,
      int detector=0
      ) {
      // logE: log10(Energy [eV])
      // Xmax: depth of maximum of em cascade [g/cm2]
      // Nmu: local relative muon number (1=proton QGSJetII)
      // theta: EAS zenith angle [rad]
      // z: altitude above sea level [m]
      // r: distance to the core [m]
      // psi: asymmetry angle in the shower plane [rad]
      // component: 0: muon, 1: em, 2: muon-induced em, 3: hadronic-induced em, -1: sum of them
      // detector: 0: WCD, 1: 2m2 ASCII, 2: 10m2 AMIGA
      return Param[detector]->GetSignal(r*100,Xmax,logE,theta,psi,rhoGround,z*100,Nmu,component,atmos);
    };
    static CDASUniv *Instance() {
      if (!_instance) _instance = new CDASUniv();
      return _instance;
    };
    UnivParamNS::UnivParam *Param[3];
  private:
    static CDASUniv * _instance;
    CDASUniv() {
      for (int i=0;i<3;i++) Param[i]=new UnivParamNS::UnivParam(i);
      atmos=-1;
      rhoGround=1.04e-3;
    };
    CDASUniv(CDASUniv const&);              // Don't Implement
    void operator=(CDASUniv const&); // Don't implement
};

CDASUniv& UnivCDAS() {
  return *CDASUniv::Instance();
}

#if 0
class UnivShower {
  private:
    double sig[4],m[4],s[4],r,azimuth,z;
    int tm;
    UnivShower() {
      InitTimeModel();
    }
  public:
    double Nmu;
    double Xmax;
    double logE;
    double Theta;
    double Phi;
    double XCore;
    double YCore;
    double ZCore;
    double GetSignal(int comp=-1) {
      return UnivCDAS().GetSignal(logE,Xmax,Nmu,Theta,z+kSTC::A0,r,azimuth,comp);
    }
    void SetStation(double x, double y, double z) {  // compute sig,m,s,r,azimuth
      double dx=x-XCore;
      double dy=y-YCore;
      double dz=z-ZCore;
      double azimuth=atan2(cos(event->fRec.Phi())*y-sin(event->fRec.Phi())*x,cos(event->fRec.Theta())*(sin(event->fRec.Phi())*y+cos(event->fRec.Phi())*x));
      for (int k=0;k<4;k++)
        sig[k]=GetSignal(k);
      if (sig[0]+sig[1]+sig[2]+sig[3]>30 && r>200 && r<2200)
        for (int k=0;k<4;k++)
          TimeModel(r,z+kSTC::A0,Xmax,Theta,azimuth,k,&m[k],&s[k]);
    }
    double GetTracePerNs(int t) {
      double ret=0;
      for (int l=0;l<4;l++)
        ret+=sig[l]*ROOT::Math::lognormal_pdf(t,m[l], s[l], 0);
      return l;
    }
};
#endif

#endif