Tools/InclinedShowers/TankResponse/USCInter/TankResponse.cc

Go to the documentation of this file.

#include "TankResponse.h"

#include <utl/AugerException.h>
#include <utl/Branch.h>
#include <utl/AugerUnits.h>
#include <utl/NormalDistribution.h>
#include <utl/Math.h>

#include <string>
#include <sstream>
#include <fstream>

#include <TFile.h>
#include <TH1D.h>

using namespace USCInterTankResponseNS;
using namespace utl;
using namespace std;
using namespace tls;

inline
double
TriggerPDF(const double signal, const double* pars)
{
  return pow(signal,pars[0])/(pow(signal,pars[0]) + pars[1]);
}

TankResponse& TankResponse::GetInstance(const Branch topBranch)
{
  static TankResponse singleton(topBranch);
  return singleton;
}

TankResponse::TankResponse(const Branch branch)
{
  branch.GetChild("Simulation").GetData(fSimulationMode);

  string dataFilePath;
  branch.GetChild("DataFile").GetData(dataFilePath);

  string dataSilentFilePath;
  branch.GetChild("DataSilentFile").GetData(dataSilentFilePath);

  string dataSilentSimFilePath;
  branch.GetChild("DataSilentSimFile").GetData(dataSilentSimFilePath);

  vector<double> buffer;
  branch.GetChild("TriggerParameters").GetData(buffer);
  fTriggerParameters[0][0] = buffer[0];
  fTriggerParameters[0][1] = buffer[1];

  vector<double>().swap(buffer);
  branch.GetChild("TriggerSimParameters").GetData(buffer);
  fTriggerParameters[1][0] = buffer[0];
  fTriggerParameters[1][1] = buffer[1];

  TFile* f = TFile::Open(dataFilePath.c_str(),"READ");
  if (!f) {
    ostringstream msg;
    msg << "ERROR: could not open file " << dataFilePath;
    throw NoDataForModelException(msg.str());
  }

  const int Nz   = 15;
  const int Nr   = 8;
  const int Nmu  = 8;
  boost::multi_array<double,3> data(boost::extents[NBIN][Nz][Nr]);
  boost::multi_array<double,2> data_avg(boost::extents[Nz][Nr]);
  boost::multi_array<double,2> data_sigma(boost::extents[Nz][Nr]);

  for (int nmu = 1; nmu <= Nmu; ++nmu)
  {
    for (int ze = 0; ze < Nz; ++ze)
    {
      int theta = 60 + 2 * ze;
      for (int rr = 0; rr < Nr; ++rr)
      {
        //open the histogram
        ostringstream histoID;
        if (nmu==1)
          histoID << "thPmt_" << theta << "_" << rr << "_0";
        else
          histoID << "thConv_" << theta << "_" << rr << "_0" << "_" << nmu;

        TH1D *histResponse = 0;
        histResponse = static_cast<TH1D*>(f->Get(histoID.str().c_str()));
        if (!histResponse)
        {
          ostringstream msg;
          msg << "ERROR: Histogram " << histoID.str() << " not found in " << dataFilePath << "\n";
          throw NoDataForModelException(msg.str());
        }

        if (nmu==1)
        {
          data_avg[ze][rr]   = histResponse->GetMean();
          data_sigma[ze][rr] = histResponse->GetRMS();
          fDS = histResponse->GetBinWidth(1);
        }

        for (int bin = 0; bin < histResponse->GetNbinsX(); bin++)
          data[bin][ze][rr] = histResponse->GetBinContent(bin+1);
      }
    }

    fUSCTankModel[nmu-1] = USCTankModel(GetThetaMin(),GetThetaMax(),GetLogrMin(),GetLogrMax(),data);

    if (nmu==1)
    {
      fUSCAvgTankModel   = USCAvgTankModel(GetThetaMin(),GetThetaMax(),GetLogrMin(),GetLogrMax(),data_avg);
      fUSCSigmaTankModel = USCSigmaTankModel(GetThetaMin(),GetThetaMax(),GetLogrMin(),GetLogrMax(),data_sigma);
    }
  }

  delete f;

  // Load silent station probability
  ifstream fOutProb;
  int dummy0, dummy1;

  fOutProb.open(dataSilentSimFilePath.c_str());
  for(int pp=0;pp<16;pp++)
    for(int k=0;k<30;k++)
      fOutProb >> dummy0 >> dummy1 >> fProbNoTriggerSilent[1][pp][k];
  fOutProb.close();

  fOutProb.open(dataSilentFilePath.c_str());
  for(int pp=0;pp<16;pp++)
    for(int k=0;k<30;k++)
      fOutProb >> dummy0 >> dummy1 >> fProbNoTriggerSilent[0][pp][k];
  fOutProb.close();
}

double
TankResponse::Mean(const double theta, const double r, const ulong nmuon)
  const
{
  double th = theta;
  if(th < GetThetaMin())
    th = GetThetaMin();

  if(th > GetThetaMax())
    th = GetThetaMax();

  double logr = log10(r/meter);

  if(logr < GetLogrMin())
    logr = GetLogrMin();

  if(logr > GetLogrMax())
    logr = GetLogrMax();

  const double a = fUSCAvgTankModel(th,logr);
  return a*nmuon;
}

double
TankResponse::StDev(const double theta, const double r, const ulong nmuon)
  const
{
  double th = theta;
  if(th < GetThetaMin())
    th = GetThetaMin();

  if(th > GetThetaMax())
    th = GetThetaMax();

  double logr = log10(r/meter);

  if(logr < GetLogrMin())
    logr = GetLogrMin();

  if(logr > GetLogrMax())
    logr = GetLogrMax();

  const double b = fUSCSigmaTankModel(th,logr);

  return b*sqrt(nmuon);
}

double
TankResponse::PDF(const double signal, const double theta, const double r, const ulong nmuon)
  const
{
  double th = theta;
  if(th < GetThetaMin())
    th = GetThetaMin();

  if(th > GetThetaMax())
    th = GetThetaMax();

  double logr = log10(r/meter);

  if(logr < GetLogrMin())
    logr = GetLogrMin();

  if(logr > GetLogrMax())
    logr = GetLogrMax();

  if (nmuon <= GetMaxMuonNumber()) {

    if (signal <= 0)
      return 0.0;

    static vector<double> interpolatedTrace(NBIN);
    fUSCTankModel[nmuon-1](interpolatedTrace, th, logr);

    size_t i = static_cast<size_t>(signal/fDS);
    if (i > NBIN-1)
      i = NBIN-1;

    // I don't understand this (HD)

    double prob = interpolatedTrace[i];
    prob *= TriggerPDF(signal, fTriggerParameters[int(fSimulationMode)]);
    if (prob < 0) prob = 0;

    return prob;

  } else {

    const double a = fUSCAvgTankModel(th,logr);
    const double b = fUSCSigmaTankModel(th,logr);

    const double mean  = a*nmuon;
    const double sigma = b*sqrt(nmuon);
    return NormalPDF(signal,mean,sigma);
  }
}

double
TankResponse::CDF(const double signalmin, const double theta, const double r, const ulong nmuon)
  const
{
  double th = theta;
  if(th < GetThetaMin())
    th = GetThetaMin();

  if(th > GetThetaMax())
    th = GetThetaMax();

  double logr = log10(r/meter);

  if(logr < GetLogrMin())
    logr = GetLogrMin();

  if(logr > GetLogrMax())
    logr = GetLogrMax();

  if (nmuon <= 30) {
    double theta_deg = th/degree;
    int c_bin = int( (theta_deg - 59.0)/2.0 );
    int entry = nmuon-1;

    // I don't understand this (HD)

    double prob_silent = 1.0-fProbNoTriggerSilent[int(fSimulationMode)][c_bin][entry];

    return prob_silent;

  } else {

    const double a = fUSCAvgTankModel(th,logr);
    const double b = fUSCSigmaTankModel(th,logr);

    double mean  = a*nmuon;
    double sigma = b*sqrt( (double)nmuon );
    return NormalCDF(signalmin, mean, sigma);
  }
}