UnivParam.cc

Go to the documentation of this file.

#include "UnivParam.h"
#include "Atmosphere.h"
#include "GeomAsym.h"
#include "AtmCorr.h"

using namespace UnivParamNS;


UnivParam::UnivParam(int DetectorType)
{
  if (DetectorType < 0 || DetectorType > 2) {
    printf("UnivParam: Wrong Detector Type %d \n", DetectorType);
    exit(1);
  }
  fDetectorType = DetectorType;

}

double UnivParam::DXFunc(double x, double logE,  double* parDX)
{

  double val = 0.;
  int icomp = int(parDX[nParDX - 1]);

  double Sref0 = parDX[0];
  double fS_e = parDX[1];
  double Sref = Sref0 *   pow(10., (logE - 19.) * fS_e);

  double DXmax = parDX[2] ;

  double L0 = parDX[3];
  double fL_e = parDX[4];
  double lambda = L0 * (1. + (logE - 19.) * fL_e) ;

  double DX0 = parDX[5];

  if (x < DX0) val = 0.;
  else val = Sref * pow((x - DX0) / (DXref[icomp] - DX0), (DXmax - DX0) / lambda) * exp((DXref[icomp] - x) / lambda);

  return val;
}

//
double UnivParam::RFuncPar(int icomp, int iparDX , int ipar)
{
  double val = 0.;
  if (iparDX == 0 && fDetectorType == 0) val = parSref_WCD[icomp][ipar];
  else if (iparDX == 1 && fDetectorType == 0) val = parGamma_WCD[icomp][ipar];
  else if (iparDX == 2 && fDetectorType == 0) val = parDXmax_WCD[icomp][ipar];
  else if (iparDX == 3 && fDetectorType == 0) val = parLambda_WCD[icomp][ipar];
  else if (iparDX == 4 && fDetectorType == 0) val = parfLambdaE_WCD[icomp][ipar];
  else if (iparDX == 5 && fDetectorType == 0) val = parDX0_WCD[icomp][ipar];
  else if (iparDX == 0 && fDetectorType == 1) val = parSref_Scin[icomp][ipar];
  else if (iparDX == 1 && fDetectorType == 1) val = parGamma_Scin[icomp][ipar];
  else if (iparDX == 2 && fDetectorType == 1) val = parDXmax_Scin[icomp][ipar];
  else if (iparDX == 3 && fDetectorType == 1) val = parLambda_Scin[icomp][ipar];
  else if (iparDX == 4 && fDetectorType == 1) val = parfLambdaE_Scin[icomp][ipar];
  else if (iparDX == 5 && fDetectorType == 1) val = parDX0_Scin[icomp][ipar];
  else if (iparDX == 0 && fDetectorType == 2) val = parSref_MD[icomp][ipar];
  else if (iparDX == 1 && fDetectorType == 2) val = parGamma_MD[icomp][ipar];
  else if (iparDX == 2 && fDetectorType == 2) val = parDXmax_MD[icomp][ipar];
  else if (iparDX == 3 && fDetectorType == 2) val = parLambda_MD[icomp][ipar];
  else if (iparDX == 4 && fDetectorType == 2) val = parfLambdaE_MD[icomp][ipar];
  else if (iparDX == 5 && fDetectorType == 2) val = parDX0_MD[icomp][ipar];
  return val;
}

//
double UnivParam::RFunc(double r, int icomp, int iparDX)
{
  double par[4];
  for (int ipar = 0; ipar < 4; ipar++) par[ipar] = RFuncPar(icomp, iparDX, ipar);
  return RFunc(r,  icomp,  iparDX , par);
}



double UnivParam::RFunc(double r, int icomp, int iparDX , double* par)
{
  double r_meter = r / 1.e2;

  if (iparDX == (nParDX - 1))
    return icomp;


  //-------------------------------------------------
  // Sref
  //-------------------------------------------------
  if (iparDX == 0)                     // NKG
    return par[0] * pow(r_meter / 1000., -par[1]) * pow((par[3] + r_meter) / (par[3] + 1000.), -par[2]);
  //-------------------------------------------------
  // Scaling
  //-------------------------------------------------
  else if (iparDX == 1) {
    // expo+cte
    return par[0] + exp(-r_meter / 500.) * par[1];
  }
  //-------------------------------------------------
  // DXmax
  //-------------------------------------------------

  else if (iparDX == 2) {
    if (icomp == 1) {  // gaisser-hillas
      double max = par[0];
      double rmax = par[1];
      double rlambda = par[2];
      double r0 = par[3];
      if (rmax < r0 || rlambda < 0.)
        return 1.e10;

      if (r_meter > r0)
        return max * pow((r_meter - r0) / (rmax - r0), (rmax - r0) / rlambda) * exp((rmax - r_meter) / rlambda);
      else
        return 0.;
    } else // powerlaw + cte
      return par[0] + par[1] * pow(r_meter / 1000., par[2]);
  }
  //-------------------------------------------------
  // Lambda
  //-------------------------------------------------
  else if (iparDX == 3) {
    if (icomp == 1)                                   // pol1 in r + expo
      return par[0] + r_meter / 1000.*par[1] - par[2] * exp(-r_meter / 150.);
    else if (icomp == 3)                              //pol1
      return par[0] + r_meter / 1000.*par[1];
    else
      return par[0] + par[1] * exp(par[2] * r_meter / 1000.); // expo + cte
  }
  //-------------------------------------------------
  // fLambda_E
  //-------------------------------------------------
  else if (iparDX == 4) {
    // pol1  in r
    if (icomp == 1) return par[0] + r_meter / 1000.*par[1];
    else            return par[0];
  }
  //-------------------------------------------------
  // DX0
  //-------------------------------------------------
  else if (iparDX == 5)
    return par[0];

  return -100.;
}



double UnivParam::GetS0(double r, double DX, double logE,  int icomp)
{
  //Signal in an spherical tank with shadow effects corrected and for a reference ground density of 1.04e-3 g/cm^3
  //   icomp=0,1  Muon/Em signal in [VEM]
  //   icomp=2,3  Ratio of Em Muon/Em Hadron to muonic signal
  //
  // r             [cm]
  // DX            [g/cm^2]    Muon/Em Muon-> DX_mu     Em/Em Had -> DX_edep
  // log(E)        [eV]
  // icomp=0,1,2,3 Muon/Em/Em Muon/Em had  -1 sum

  if (r > rLimit_high) r = rLimit_high;
  if (r < rLimit_low) r = rLimit_low;

  double parDX[nParDX];
  for (int iparDX = 0; iparDX < nParDX; iparDX++)
    parDX[iparDX] = RFunc(r, icomp, iparDX);

  double S0 = DXFunc(DX, logE, parDX);
  return S0;
}

double UnivParam::GetS0(double r, double DX, double logE, int icomp, double Nmu)
{
  if (icomp < 0 || icomp > 3) return 0.;
  if (r > rLimit_high) r = rLimit_high;
  if (r < rLimit_low) r = rLimit_low;

  double val = GetS0(r, DX, logE, icomp);

  double alpha = GetAlphaFluct(r, icomp);
  if (alpha == UNDEF) return UNDEF;
  return val * (1. + (Nmu - 1.) * alpha);
}


double UnivParam::GetSignal(double r,  double Xmax_Edep , double logE, double theta, double psi,  double rhoGround, double hGround, double Nmu,  int icomp0, int iatm)
{

  if (isnan(Xmax_Edep) || Xmax_Edep < 0.) return 0.;
  AtmosphereNS::Atmosphere fatm;
  if (iatm > 0 && iatm < 13) fatm.SetCurrentAtmosphere(iatm);

  if (hGround < 0.)   hGround = AtmosphereNS::hGroundRef;
  if (rhoGround < 0.)   rhoGround = fatm.GetDensity(hGround);

  double DX[4];
  for (int icomp = 0; icomp < 4; icomp++)
    DX[icomp] = fatm.Get_DX_DL(r,  psi,  Xmax_Edep + XmaxShift[icomp]   , theta,  hGround, UseDL[icomp], UseDiffusive[icomp]);

  return GetSignal(r, DX, logE, theta, psi, rhoGround, hGround, Nmu, icomp0);

}

double UnivParam::GetSignal(double r,  double* DX , double logE, double theta, double psi,  double rhoGround, double hGround, double Nmu,  int icomp0)
{
  // r             [cm]
  // Xmax_edep     [g/cm²]
  // log(E)        [eV]
  // theta/psi     [rad]
  // rhoGround     [g/cm^3]
  // hGround       [cm]
  // Nmu           relative local Nmu
  // icomp=0,1,2,3 Muon/Em/Em Muon/Em had  -1 sum

  if (isnan(r) || isnan(logE) || isnan(theta)  || isnan(Nmu)) return 0;
  if (icomp0 < -1 || icomp0 > 3) return 0.;
  if (fDetectorType == 2 && !(icomp0 == 0 || icomp0 == -1))    return 0.;
  if (hGround <= 0.)                                      return 0.;
  if (rhoGround <= 0.)                                      return 0.;

  if (Nmu < 0.) Nmu = 0.;
  if (Nmu > 3.5) Nmu = 3.5;
  if (r > rLimit_high) r = rLimit_high;
  if (r < rLimit_low) r = rLimit_low;
  if (theta > 65 * M_PI / 180.)               theta = 65.*M_PI / 180.;


  AtmCorrNS::AtmCorr fatmcorr(fDetectorType);
  GeomAsymNS::GeomAsym fgeom(fDetectorType);

  double sum = 0.;
  for (int icomp = 0; icomp < 4; icomp++) {
    if (!(icomp == icomp0 || icomp0 == -1)) continue;
    if (fDetectorType == 2 && icomp > 0) continue;

    double DX0 = RFunc(r, icomp, 5);
    if (DX[icomp] <= DX0)               continue;
    if (!UseDiffusive[icomp] && DX[icomp] <= 0.) continue;

    double S0 = GetS0(r,  DX[icomp],   logE,   icomp , Nmu)  ;
    double S =  S0 *  fgeom.GetGeomAsym(DX[icomp],    r,  theta,  psi,  icomp)
                /  fatmcorr.GetCorrTypeI(r, icomp, rhoGround, hGround) /  fatmcorr.GetCorrTypeII(r, theta, icomp, hGround) / fatmcorr.GetCorrTypeIII(r, theta, psi, icomp);

    if (S <= 0.) S = 0.;
    sum += S;
  }
  return sum;
}


double UnivParam::GetLogE(double vem, double r, double Xmax_Edep , double theta, double psi, double rhoGround, double hGround, double Nmu, int iatm)
{
  double logE = 17.00;
  double vemG = GetSignal(r, Xmax_Edep, logE, theta, psi, rhoGround, hGround, Nmu, -1, iatm);
  if (vem < vemG) return logE;

  double dlogE = 0.05;
  while (fabs(vem / vemG - 1.) > 0.0001) {
    if (vem > vemG) logE += dlogE;
    else {
      logE -= dlogE;
      dlogE /= 2.;
    }
    vemG = GetSignal(r, Xmax_Edep, logE, theta, psi, rhoGround, hGround, Nmu, -1, iatm);
  }
  return logE;
}

double UnivParam::GetAlphaFluct(double r, int icomp)
{
  if (icomp == 0 || icomp == 2) return 1.;
  else if (icomp == 1) return -0.075;
  else if (icomp == 3) {
    double r_meter = r / 1.e2;
    double alpha0 = 1.12, alpha1 = 1.25, B = 6.;

    double C = (alpha1 - alpha0) / (exp(-B) - 1.0);
    double A = alpha0 - C;

    return A + C * exp(-B * r_meter / 1.e3);
  }
  return UNDEF;

}