AachenTable.cc

Go to the documentation of this file.

#include "AachenTable.h"

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

#include <TFile.h>
#include <TH1D.h>
#include <TF1.h>
#include <TGraph.h>
#include <TCanvas.h>

#include <iostream>
#include <sstream>
#include <algorithm>

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

AachenTable::AachenTable()
{}

AachenTable::AachenTable(const string& fileName)
{
  OpenFile(fileName);
}

AachenTable::~AachenTable()
{
}

void AachenTable::OpenFile(const string& fileName)
{
  TFile* f = TFile::Open(fileName.c_str(),"READ");
  if (!f) {
    ostringstream msg;
    msg << "ERROR: could not open file " << fileName;
    throw NoDataForModelException(msg.str());
  }

  cout << "AachenTankResponse convolving";
  for (int theta = GetThetaMinDeg(); theta <= GetThetaMaxDeg(); theta+=2) {
    cout << "." << flush;

    fTable.push_back(TableItem());
    TableItem& item = fTable.back();

    item.theta = theta*degree;

    //read the histogram
    ostringstream histoID;
    histoID << "p" << theta << "_1";

    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 " << fileName << "\n";
      throw NoDataForModelException(msg.str());
    }

    // Rough fix to reduce the failure rate of MINOS,
    // which gets confused by zero probabilities in rare cases
    double firstMean = histResponse->GetMean();
    double firstSigma = histResponse->GetRMS();
    for (int ix=1;ix<histResponse->GetNbinsX()+1;++ix)
    {
      double y = histResponse->GetBinContent(ix);
      if (y==0) {
        const double x = histResponse->GetXaxis()->GetBinCenter(ix);
        histResponse->SetBinContent(ix,NormalPDF(x,firstMean,firstSigma));
      }
    }

    // fill the first TPF
    item.GetTPF(1).SetScale(
      histResponse->GetXaxis()->GetXmin()+histResponse->GetXaxis()->GetBinWidth(1)/2,
      histResponse->GetXaxis()->GetBinWidth(1));
    for (int ix=1;ix<histResponse->GetNbinsX()+1;++ix)
    {
      const double& y = histResponse->GetBinContent(ix);
      item.GetTPF(1).PushBack(y);
    }
    item.GetTPF(1).Normalise();
    item.GetMean(1) = item.GetTPF(1).Moment(1);

    // auto-convolution
    TH1D* h0 = histResponse;
    const int n = h0->GetNbinsX();
    double* conv = new double[n];
    double* last = new double[n];
    for (int i=0;i<n;++i) last[i] = h0->GetBinContent(i+1);

    for (int nmu=2;nmu<GetTableLimit();++nmu) {
      for (int i=0;i<n;++i) {
        double psum = 0;
        for (int ii=0;ii<=i;++ii) {
          const double p1 = last[ii];
          const double p2 = h0->GetBinContent(i-ii+1);
          psum += p1*p2;
        }
        conv[i] = psum;
      }

      // fill the next TPF
      item.GetTPF(nmu).SetScale(
        histResponse->GetXaxis()->GetXmin()+histResponse->GetXaxis()->GetBinWidth(1)/2,
        histResponse->GetXaxis()->GetBinWidth(1));
      for (int i = 0;i<n;++i)
      {
        const double& y = conv[i];
        item.GetTPF(nmu).PushBack(y);
      }
      item.GetTPF(nmu).Normalise();
      item.GetMean(nmu) = item.GetTPF(nmu).Median();

      for (int i=0;i<n;++i) last[i] = conv[i];
    }

    delete histResponse;

    // calculate extrapolation parameters
    TGraph gMean;
    TGraph gSigma;
    for (int nmu=1;nmu<GetTableLimit();++nmu) {
      const double mean = item.GetTPF(nmu).Moment(1);
      /* use truncated variance to get better control over the long tails */
      const double truncVar = item.GetTPF(nmu).Variance(item.GetTPF(nmu).Quantile(0.05),item.GetTPF(nmu).Quantile(0.95));
      gMean.SetPoint(gMean.GetN(),nmu,mean);
      gSigma.SetPoint(gSigma.GetN(),nmu,sqrt(truncVar));
    }

    TF1 fitMean("fitMean","[0]*x");
    gMean.Fit(&fitMean,"0Q");
    item.parMean = fitMean.GetParameter(0);

    TF1 fitSigma("fitSigma","sqrt([0]+[1]*x)");
    gSigma.Fit(&fitSigma,"0Q");
    item.parSigma[0] = fitSigma.GetParameter(0);
    item.parSigma[1] = fitSigma.GetParameter(1);
  }
  cout << endl;

  delete f;
}

double AachenTable::PDF(double signal, double theta, int nmu) const
{
  if (theta <= fTable[0].theta) {
    return fTable[0].GetTPF(nmu).Density(signal);
  }

  if (theta >= fTable.back().theta) {
    return fTable.back().GetTPF(nmu).Density(signal);
  }

  const InternalTable::const_iterator itu = UpperLimit(theta);
  const TableItem& upper = *itu;
  const TableItem& lower = *(itu-1);

  const double m = (theta - lower.theta)/(upper.theta-lower.theta);

  const double delta_mean = signal - (lower.GetMean(nmu) + m*(upper.GetMean(nmu)-lower.GetMean(nmu)));

  const double p0 = lower.GetTPF(nmu).Density(delta_mean + lower.GetMean(nmu));
  const double p1 = upper.GetTPF(nmu).Density(delta_mean + upper.GetMean(nmu));

  return p0 + m*(p1-p0);
}

double AachenTable::CDF(double signalmin, double theta, int nmu) const
{
  if (theta <= fTable[0].theta)
    return 1.0 - fTable[0].GetTPF(nmu).Integral(signalmin);

  if (theta >= fTable.back().theta)
    return 1.0 - fTable.back().GetTPF(nmu).Integral(signalmin);

  const InternalTable::const_iterator itu = UpperLimit(theta);
  const TableItem& upper = *itu;
  const TableItem& lower = *(itu-1);

  const double m = (theta - lower.theta)/(upper.theta-lower.theta);

  const double p0 = 1.0 - lower.GetTPF(nmu).Integral(signalmin);
  const double p1 = 1.0 - upper.GetTPF(nmu).Integral(signalmin);

  return p0 + m*(p1-p0);
}

double AachenTable::Moment(double theta, int nmu, int order) const {
  if (theta <= fTable[0].theta)
    return fTable[0].GetTPF(nmu).Moment((unsigned short)order);

  if (theta >= fTable.back().theta)
    return fTable.back().GetTPF(nmu).Moment((unsigned short)order);

  const InternalTable::const_iterator itu = UpperLimit(theta);
  const TableItem& upper = *itu;
  const TableItem& lower = *(itu-1);

  const double m = (theta - lower.theta)/(upper.theta-lower.theta);

  const double m0 = lower.GetTPF(nmu).Moment((unsigned short)order);
  const double m1 = upper.GetTPF(nmu).Moment((unsigned short)order);

  return m0 + m*(m1-m0);
}

double AachenTable::GetFittedMean(double theta, int nmu) const {
  if (theta <= fTable[0].theta)
    return fTable[0].GetFittedMean(nmu);

  if (theta >= fTable.back().theta)
    return fTable.back().GetFittedMean(nmu);

  const InternalTable::const_iterator itu = UpperLimit(theta);
  const TableItem& upper = *itu;
  const TableItem& lower = *(itu-1);

  const double m = (theta - lower.theta)/(upper.theta-lower.theta);

  const double m0 = lower.GetFittedMean(nmu);
  const double m1 = upper.GetFittedMean(nmu);

  return m0 + m*(m1-m0);
}

double AachenTable::GetFittedSigma(double theta, int nmu) const {
  if (theta <= fTable[0].theta)
    return fTable[0].GetFittedSigma(nmu);

  if (theta >= fTable.back().theta)
    return fTable.back().GetFittedSigma(nmu);

  const InternalTable::const_iterator itu = UpperLimit(theta);
  const TableItem& upper = *itu;
  const TableItem& lower = *(itu-1);

  const double m = (theta - lower.theta)/(upper.theta-lower.theta);

  const double m0 = lower.GetFittedSigma(nmu);
  const double m1 = upper.GetFittedSigma(nmu);

  return m0 + m*(m1-m0);
}