FDetector/Pixel.cc

Go to the documentation of this file.

#include <boost/lexical_cast.hpp>

#include <det/VManager.h>
#include <det/Detector.h>

#include <fdet/Pixel.h>
#include <fdet/Camera.h>
#include <fdet/Eye.h>
#include <fdet/Telescope.h>
#include <fdet/FDetector.h>

#include <utl/CoordinateSystem.h>
#include <utl/AugerException.h>
#include <utl/MathConstants.h>
#include <utl/Math.h>

#include <fwk/CoordinateSystemRegistry.h>

using namespace fdet;
using namespace det;
using namespace std;
using namespace utl;
using namespace fwk;


Pixel::Pixel(const unsigned int eyeId, const unsigned int telescopeId,
             const unsigned int pixelId,
             const std::string& physEyeIdString,
             const std::string& physTelIdString) :
  //fDirection()
  //fCoordinateSystem()
  //fWavelengthCalib
  fTimeOffset(0),
  fPixelStatus(Pixel::eUnknown),
  fPixelOpticalEfficiencyStatus(Pixel::eUnknown),
  fSimulatedEndToEnd(0),
  fEyeId(eyeId),
  fTelescopeId(telescopeId),
  fPixelId(pixelId),
  fPhysicalEyeIdString(physEyeIdString),
  fPhysicalTelescopeIdString(physTelIdString),
  fSolidAngle(0),
  fQEfficiency(0)
{
  fPixelIdString = boost::lexical_cast<string>(pixelId);
}


Pixel::~Pixel()
{
  delete fTimeOffset;
  delete fSimulatedEndToEnd;
  delete fSolidAngle;
  delete fQEfficiency;
}


void
Pixel::Update()
{
  delete fTimeOffset;
  delete fSimulatedEndToEnd;
  delete fSolidAngle;
  delete fQEfficiency;
  fDirection = utl::Vector();
  fCoordinateSystem.reset();
  fTimeOffset = 0;
  // pixel status is cached via fdet::Telescope!!
  //  fPixelStatus = Pixel::eUnknown;
  fSimulatedEndToEnd = 0;
  fSolidAngle = 0;
  fQEfficiency = 0;
}


unsigned int
Pixel::GetChannelId()
  const
{
  return GetTelescope().GetCamera().GetChannelId(fPixelId);
}


const Vector&
Pixel::GetDirection()
  const
{
  if (!fDirection.GetMag2())
    fDirection = GetTelescope().GetCameraPixelDirection(GetId());
  return fDirection;
}


const CoordinateSystemPtr&
Pixel::GetPixelCoordinateSystem()
  const
{
  if (!fCoordinateSystem) {
    const Telescope& telDet = GetTelescope();
    const CoordinateSystemPtr telCS = telDet.GetTelescopeCoordinateSystem();
    const Camera& camera = telDet.GetCamera();

    const double rFocal = camera.GetRadiusFocal();
    const double heightMercedes = camera.GetMercedesHeight();
    const double pixelVertex = 26.34*mm;

    const Vector direction = -GetDirection();

    const double fat = sqrt (Sqr(rFocal-heightMercedes) - Sqr(pixelVertex));

    CoordinateSystemPtr auxCS = CoordinateSystem::Translation(direction*fat, telCS);
    CoordinateSystemPtr backCS = CoordinateSystem::RotationY(180.*deg, auxCS);

    const double phi_rot = direction.GetPhi(backCS);
    const double theta_rot = direction.GetTheta(backCS);

    CoordinateSystemPtr auxCS2 = CoordinateSystem::RotationZ(phi_rot, backCS);
    CoordinateSystemPtr auxCS3 = CoordinateSystem::RotationY(theta_rot, auxCS2);
    fCoordinateSystem = CoordinateSystem::RotationZ(kPi/2 - phi_rot, auxCS3);
  }
  return fCoordinateSystem;
}


double
Pixel::GetSimulatedDiaPhoton2ADC(const string& configSignature, const double wavelength)
  const
{
  const Telescope& tel = GetTelescope();
  const CoordinateSystemPtr& telCS = tel.GetTelescopeCoordinateSystem();
  const double cosTheta = GetDirection().GetCosTheta(telCS);
  const double calib = GetSimulatedEndToEndCalibration(configSignature);
  return tel.GetModelRelativeEfficiency(wavelength) / (calib * cosTheta);
}


double
Pixel::GetEndToEndCalibrationAtReferenceWavelength()
  const
{
  const double refWl = Detector::GetInstance().GetFDetector().GetReferenceLambda();
  return GetEndToEndCalibration(refWl);
}


double
Pixel::GetEndToEndCalibrationConstantAtReferenceWavelength()
  const
{
  const double refWl = Detector::GetInstance().GetFDetector().GetReferenceLambda();
  return GetEndToEndCalibrationConstant(refWl);
}


double
Pixel::GetOpticalEfficiencyCorrectionAtReferenceWavelength()
  const
{
  const double refWl = Detector::GetInstance().GetFDetector().GetReferenceLambda();
  return GetOpticalEfficiencyCorrection(refWl);
}


const TabulatedFunction&
Pixel::GetEndToEndCalibrationConstant()
  const
{
  // Pass job of filling all calib consts to the parent telescopes
  GetTelescope().CachePixelCalibrations();
  return fWavelengthCalib;
}


const TabulatedFunction&
Pixel::GetOpticalEfficiencyCorrection()
  const
{
  // Pass job of filling all optical efficiency corrections to the parent telescopes
  GetTelescope().CachePixelOpticalEfficiencyCorrections();
  return fWavelengthOpticalEfficiencyCorrection;
}


//for particular wavelength we multiply calibration constant by optical efficiency correction
double
Pixel::GetEndToEndCalibration(const double wavelength)
  const
{
  const Telescope& telDet = GetTelescope();
  telDet.CachePixelCalibrations();
  telDet.CachePixelOpticalEfficiencyCorrections();

  double opticalEfficiencyCorrection = 0;
  if (wavelength < fWavelengthOpticalEfficiencyCorrection.XFront()) {
    opticalEfficiencyCorrection = fWavelengthOpticalEfficiencyCorrection.YFront();
  } else if (wavelength > fWavelengthOpticalEfficiencyCorrection.XBack()) {
    opticalEfficiencyCorrection = fWavelengthOpticalEfficiencyCorrection.YBack();
  } else {
    opticalEfficiencyCorrection = fWavelengthOpticalEfficiencyCorrection(wavelength);
  }

  const double wavelengthFront = fWavelengthCalib.XFront();
  if (wavelength < wavelengthFront) {

    // wavelength outside (below) of the absolute calibration boundaries
    // => extrapolate with the relative calibration
    return fWavelengthCalib.YFront() / telDet.GetMeasuredRelativeEfficiency(wavelengthFront) * opticalEfficiencyCorrection;
  }

  const double wavelengthBack = fWavelengthCalib.XBack();
  if (wavelength > wavelengthBack) {

    // ... above ...
    return fWavelengthCalib.YBack() / telDet.GetMeasuredRelativeEfficiency(wavelengthBack) * opticalEfficiencyCorrection;
  }

  // wavelength is within tabulated function boundaries => interpolate
  return fWavelengthCalib(wavelength) * opticalEfficiencyCorrection;
}


double
Pixel::GetEndToEndCalibrationConstant(const double wavelength)
  const
{
  const Telescope& telDet = GetTelescope();
  telDet.CachePixelCalibrations();

  const double wavelengthFront = fWavelengthCalib.XFront();
  if (wavelength < wavelengthFront) {

    // wavelength outside (below) of the absolute calibration boundaries
    // => extrapolate with the relative calibration
    return fWavelengthCalib.YFront() / telDet.GetMeasuredRelativeEfficiency(wavelengthFront);
  }

  const double wavelengthBack = fWavelengthCalib.XBack();
  if (wavelength > wavelengthBack) {

    // ... above ...
    return fWavelengthCalib.YBack() / telDet.GetMeasuredRelativeEfficiency(wavelengthBack);
  }

  // wavelength is within tabulated function boundaries => interpolate
  return fWavelengthCalib(wavelength);
}


double
Pixel::GetOpticalEfficiencyCorrection(const double wavelength)
  const
{
  const Telescope& telDet = GetTelescope();
  telDet.CachePixelOpticalEfficiencyCorrections();

  if (wavelength < fWavelengthOpticalEfficiencyCorrection.XFront()) {
    return fWavelengthOpticalEfficiencyCorrection.YFront();
  } else if (wavelength > fWavelengthOpticalEfficiencyCorrection.XBack()) {
    return fWavelengthOpticalEfficiencyCorrection.YBack();
  }

  // wavelength is within tabulated function boundaries => interpolate
  return fWavelengthOpticalEfficiencyCorrection(wavelength);
}


double
Pixel::GetDiaPhoton2PEFactor(const double wavelength, const std::string& configSignature)
  const
{

  const Telescope& telDet = GetTelescope();
  const Camera& camera = telDet.GetCamera();
  double eff = 0.;
  if (configSignature.empty())
    eff = GetSimulatedDiaPhoton2ADC(telDet.GetConfigSignature("TelescopeSimulatorKG"),
                                    wavelength);
  else 
    eff = GetSimulatedDiaPhoton2ADC(configSignature, wavelength);
  
  return eff / camera.GetElectronicsGain();
  /*
  const Telescope& tel = GetTelescope();
  const CoordinateSystemPtr& telCS = tel.GetTelescopeCoordinateSystem();
  const double cosTheta = GetDirection().GetCosTheta(telCS);
  const double calib = GetEndToEndCalibration(wavelength);
  const double diaPhoton2ADC = 1 / (calib * cosTheta);
  const double diaPhoton2PE =  diaPhoton2ADC / tel.GetCamera().GetElectronicsGain();
  return diaPhoton2PE;
  */
}


double
Pixel::GetSimulatedEndToEndCalibration(const string& configSignature)
  const
{
  return GetPixelData(fSimulatedEndToEnd, "sim_pixel_calibs", "fd_sim_calib",
                      "simulated pixel calibration", configSignature);
}


/*
double
Pixel::GetSimulatedEndToEndCalibration(const std::string& configSignature, const double wl)
  const
{
  return GetSimulatedEndToEndCalibration(configSignature);
  return fWavelengthCalib.YFront() / telDet.GetMeasuredRelativeEfficiency(wavelengthFront);
}
*/


Pixel::Status
Pixel::GetStatus()
  const
{
  // Pass job of filling pixel status to the parent telescopes
  GetTelescope().CachePixelCalibrations();
  return fPixelStatus;
}

Pixel::Status
Pixel::GetOpticalEfficiencyStatus()
  const
{
  // Pass job of filling pixel optical efficiency status to the parent telescopes
  GetTelescope().CachePixelOpticalEfficiencyCorrections();
  return fPixelOpticalEfficiencyStatus;
}


double
Pixel::GetTimeOffset()
  const
{
  return GetPixelData(fTimeOffset, "pixel_t_offset", "fd_calib", "pixel time offset");
}


unsigned int
Pixel::GetRow(const unsigned int pixelId)
  const
{
  const Telescope& tel = GetTelescope();
  return (pixelId - tel.GetFirstPixelId()) % tel.GetLastRow() + tel.GetFirstRow();
}


unsigned int
Pixel::GetColumn(const unsigned int pixelId)
  const
{
  const Telescope& tel = GetTelescope();
  return (pixelId - tel.GetFirstPixelId()) / tel.GetLastRow() + tel.GetFirstColumn();
}


const TabulatedFunction&
Pixel::GetQEfficiency()
  const
{
  return GetPixelData(fQEfficiency, "qEfficiency", "pixel", "PMT quantum efficiciency");
}


double
Pixel::GetSolidAngle()
  const
{
  if (!fSolidAngle) {
    vector<double> tmpSolidAngleList;
    GetPixelData(tmpSolidAngleList, "pixelSolidAngle", "camera", "PMT solid angle");
    fSolidAngle = new double(tmpSolidAngleList[fPixelId - 1]);
  }
  return *fSolidAngle;
}


int
Pixel::GetCloudIndex()
  const
{
  if (!fCloudValidityStamp.IsValid()) {
    std::ostringstream err;
    err << "no cloud index present (try with HasCloudFraction())";
    throw utl::DataNotFoundInDBException(err.str());
  }

  return fCloudIndex;
}


float
Pixel::GetCloudFraction()
  const
{
  const float indexToFraction[6] = {
    0.00, // 0-10
    0.20, // 10-30
    0.40, // 30-50
    0.60, // 50-70
    0.80, // 70-90
    1.00, // 90-100
  };

  return indexToFraction[GetCloudIndex()];
}


bool
Pixel::HasCloudFraction()
  const
{
  if (fCloudValidityStamp.IsValid())
    return true;
  else
    return GetTelescope().CachePixelCloudData();
}


const Telescope&
Pixel::GetTelescope()
  const
{
  const Telescope& tel =
    Detector::GetInstance().GetFDetector().GetEye(fEyeId).GetTelescope(fTelescopeId);
  return tel;
}


template<typename T>
const T& Pixel::GetPixelData(T*& requestedData,
                             const std::string& property,
                             const std::string& component,
                             const std::string& errorMsg,
                             const std::string& extraIndex)
  const
{
  if (!requestedData) {
    requestedData = new T;
    GetPixelData(*requestedData, property, component, errorMsg, extraIndex);
  }
  return *requestedData;
}


template<typename T>
void Pixel::GetPixelData(T& requestedData,
                         const std::string& property,
                         const std::string& component,
                         const std::string& errorMsg,
                         const std::string& extraIndex)
  const
{
  const det::VManager& manager = det::Detector::GetInstance().GetFManagerRegister();
  det::VManager::IndexMap indexMap;
  indexMap["pixelId"] = fPixelIdString;
  indexMap["eyeId"] = fPhysicalEyeIdString;
  indexMap["telescopeId"] = fPhysicalTelescopeIdString;
  if (!extraIndex.empty())
    indexMap["signature"] = extraIndex;

  const det::VManager::Status foundFlag =
    manager.GetData(requestedData, property, component, indexMap);

  // ERROR
  if (foundFlag == det::VManager::eNotFound) {
    std::ostringstream err;
    err << "Did not find requested component for " << errorMsg<< "; "
        << det::VManager::QueryInfoMessage(property, component, indexMap);
    ERROR(err);
    throw utl::NonExistentComponentException(err.str());
  }
}


template<typename T>
const T&
Pixel::GetPixelDataDiaphragm(T*& requestedData,
                             const std::string& property,
                             const std::string& component,
                             const std::string& errorMsg)
  const
{
  if (!requestedData) {
    requestedData = new T;
    GetPixelDataDiaphragm(*requestedData, property, component, errorMsg);
  }
  return *requestedData;
}


template<typename T>
void
Pixel::GetPixelDataDiaphragm(T& requestedData,
                             const std::string& property,
                             const std::string& component,
                             const std::string& errorMsg)
  const
{
  const det::VManager& manager = det::Detector::GetInstance().GetFManagerRegister();
  det::VManager::IndexMap indexMap;
  indexMap["pixelId"] = fPixelIdString;
  indexMap["eyeId"] = fPhysicalEyeIdString;
  indexMap["telescopeId"] = fPhysicalTelescopeIdString;
  const fdet::FDetector& detFD = Detector::GetInstance().GetFDetector();
  const bool hasCorrectorRing =
    detFD.GetEye(fEyeId).GetTelescope(fTelescopeId).HasCorrectorRing();
  indexMap["correctorId"] = (hasCorrectorRing ? "large" : "small");

  const det::VManager::Status foundFlag =
    manager.GetData(requestedData, property, component, indexMap);

  // ERROR
  if (foundFlag == det::VManager::eNotFound) {

    std::ostringstream err;
    err << "Did not find requested component for " << errorMsg<< "; "
        << det::VManager::QueryInfoMessage(property, component, indexMap);
    ERROR(err);
    throw utl::NonExistentComponentException(err.str());
  }
}


// Configure (x)emacs for this file ...
// Local Variables:
// mode: c++
// compile-command: "make -C .. FDetector/Pixel.o -k"
// End: