FParametricOpticalEfficiencyLossManager.cc

Go to the documentation of this file.

#include <string>
#include <sstream>
#include <set>
#include <map>

#include <boost/lexical_cast.hpp>

#include <fdet/FParametricOpticalEfficiencyLossManager.h>
#include <fwk/CentralConfig.h>
#include <det/Detector.h>

#include <utl/ErrorLogger.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/Reader.h>
#include <utl/TabulatedFunction.h>
#include <utl/AugerException.h>
#include <utl/ShadowPtr.h>

#include <fdet/FManagerRegister.h>
#include <fdet/FDetector.h>
#include <fdet/Telescope.h>

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


REGISTER_F_MANAGER("FParametricOpticalEfficiencyLossManager", FParametricOpticalEfficiencyLossManager);


void
FParametricOpticalEfficiencyLossManager::Init(const std::string& configLink)
{
  fHasDefaultOpticalEfficiencyCorrection = false;
  VManager::Init(configLink);
  Initialize();
  ostringstream info;
  if (fHasDefaultOpticalEfficiencyCorrection)
    info << "fixed optical efficiency correction is " << fOpticalEfficiencyCorrection;
  else
    info << " no default fixed optical efficiency correction";

  if (!fPerTelescopeOpticalEfficiencyCorrections.empty()) {
    for (map<int, map<int, vector<double> > >::const_iterator eyeIt = fPerTelescopeOpticalEfficiencyCorrections.begin(),
         eyeEnd = fPerTelescopeOpticalEfficiencyCorrections.end(); eyeIt != eyeEnd; ++eyeIt) {

      for (map<int, vector<double> >::const_iterator telIt = eyeIt->second.begin(),
           telEnd = eyeIt->second.end(); telIt != telEnd; ++telIt) {

        info << "\n    Override eye " << eyeIt->first << " tel " << telIt->first
             << ": \n epochs :";
        for (vector<unsigned int>::const_iterator epochIt = fPerTelescopeOpticalEfficiencyEpochs[eyeIt->first][telIt->first].begin(),
             epochEnd = fPerTelescopeOpticalEfficiencyEpochs[eyeIt->first][telIt->first].end(); epochIt != epochEnd; ++epochIt) {
          info << " " << *epochIt;
        }
        info << "\n corrections :";
        for (vector<double>::const_iterator corrIt = telIt->second.begin(),
             corrEnd = telIt->second.end(); corrIt != corrEnd; ++corrIt) {
          info << " " << *corrIt;
        }

      } // end for tels
    } // end for eyes
  } // end if have per-tel optical efficiency correction

  INFO(info);
}


void
FParametricOpticalEfficiencyLossManager::Initialize()
{
  Branch defaultOpticalEffBranch = fBranch.GetChild("fixedOpticalEfficiencyCorrection");
  if (defaultOpticalEffBranch) {
    defaultOpticalEffBranch.GetData(fOpticalEfficiencyCorrection);
    fHasDefaultOpticalEfficiencyCorrection = true;
  } else
    fHasDefaultOpticalEfficiencyCorrection = false;

  const Branch& perTel = fBranch.GetChild("perTelescope");
  if (perTel) {

    Branch opticalEffBranch = perTel.GetFirstChild();
    while (opticalEffBranch) {

      const AttributeMap& attrs = opticalEffBranch.GetAttributes();
      AttributeMap::const_iterator attrEnd = attrs.end();
      AttributeMap::const_iterator eyeIdIt = attrs.find("eyeId");
      AttributeMap::const_iterator telIdIt = attrs.find("telId");
      if (eyeIdIt == attrEnd || telIdIt == attrEnd) {
        ERROR("telescopeOpticalEfficiencyCorrection tag is missing eyeId or telId attributes!");
        break;
      }
      const int eyeId = boost::lexical_cast<int>(eyeIdIt->second);
      const int telId = boost::lexical_cast<int>(telIdIt->second);

      vector<unsigned int> epochs;
      vector<double> corrections;
      opticalEffBranch.GetChild("opticalEfficiencyEpochs").GetData(epochs);
      opticalEffBranch.GetChild("opticalEfficiencyCorrections").GetData(corrections);
      fPerTelescopeOpticalEfficiencyEpochs[eyeId][telId] = epochs;
      fPerTelescopeOpticalEfficiencyCorrections[eyeId][telId] = corrections;

      opticalEffBranch = opticalEffBranch.GetNextSibling();

    }

  } // end have optical eff correction per-tel

  // list of properties
  //fAvailableComponents.insert("pixel_t_offset");
  fAvailableComponents.insert("pixel_optical_efficiency_corrections");
  fAvailableComponents.insert("status");
  fAvailableComponents.insert("start_time");
  fAvailableComponents.insert("end_time");
}


VManager::Status
FParametricOpticalEfficiencyLossManager::GetData(double& /*returnData*/,
                                                 const string& /*componentProperty*/,
                                                 const string& /*componentName*/,
                                                 const IndexMap& )
  const
{
  // this manager only knows about fd_optical_efficiency
  /*if (componentName != "fd_optical_efficiency")
    return eNotFound;

  // check property availability
  if (fAvailableComponents.find(componentProperty) == fAvailableComponents.end())
    return eNotFound;
  */
  return eNotFound;
}


VManager::Status
FParametricOpticalEfficiencyLossManager::GetData(map<int, utl::TabulatedFunction>& returnData,
                                                 const string& componentProperty,
                                                 const string& componentName,
                                                 const IndexMap& componentIndex)
  const
{
  // this manager only knows about fd_optical_efficiency
  if (componentName != "fd_optical_efficiency")
    return eNotFound;

  // check property availability
  if (fAvailableComponents.find(componentProperty) == fAvailableComponents.end())
    return eNotFound;

  returnData.clear();

  if (componentProperty == "pixel_optical_efficiency_corrections") {
    const fdet::FDetector& fdet = det::Detector::GetInstance().GetFDetector();
    const double wavelength = fdet.GetReferenceLambda();

    const string& eyeIdString = componentIndex.find("eyeId")->second;
    const string& telescopeIdString = componentIndex.find("telescopeId")->second;
    const int eyeId = boost::lexical_cast<int>(eyeIdString);
    const int telId = boost::lexical_cast<int>(telescopeIdString);

    double opticalEff = 0;
    if (GetOpticalEfficiencyCorrectionForTelescope(eyeId, telId, opticalEff)) {

      for (int iPix = 1; iPix <= 440; ++iPix) {
        TabulatedFunction effVsWavelength;
        effVsWavelength.PushBack(wavelength, opticalEff);
        returnData.insert(make_pair(iPix, effVsWavelength));
      }

      return eFound;
    }
  }

  return eNotFound;
}


bool
FParametricOpticalEfficiencyLossManager::GetOpticalEfficiencyCorrectionForTelescope(const unsigned int eyeId,
                                                                                    const unsigned int telId,
                                                                                    double& value)
  const
{
  map<int, map<int, vector<double> > >::const_iterator eyeIt = fPerTelescopeOpticalEfficiencyCorrections.find(eyeId);
  if (eyeIt == fPerTelescopeOpticalEfficiencyCorrections.end()) {
    if (fHasDefaultOpticalEfficiencyCorrection) {
      value = fOpticalEfficiencyCorrection;
      return true;
    } else
      return false;
  }

  map<int, vector<double> >::const_iterator telIt = eyeIt->second.find(telId);
  if (telIt == eyeIt->second.end()) {
    if (fHasDefaultOpticalEfficiencyCorrection) {
      value = fOpticalEfficiencyCorrection;
      return true;
    } else
      return false;
  }

  map<int, map<int, vector<unsigned int> > >::const_iterator eyeEpochIt = fPerTelescopeOpticalEfficiencyEpochs.find(eyeId);
  map<int, vector<unsigned int> >::const_iterator telEpochIt = eyeEpochIt->second.find(telId);

  const TimeStamp detTime = Detector::GetInstance().GetTime();
  const unsigned int gpsSecond = detTime.GetGPSSecond();
  int thisBin = -1;
  //for (unsigned int i = 0; i < fPerTelescopeOpticalEfficiencyEpochs[eyeId][telId].size(); ++i) {
  for (unsigned int i = 0; i < telEpochIt->second.size(); ++i) {
    //if (fPerTelescopeOpticalEfficiencyEpochs[eyeId][telId][i] > gpsSecond) {
    if (telEpochIt->second[i] > gpsSecond) {
      thisBin = i - 1;
      break;
    }
  }

  //int lastBin = fPerTelescopeOpticalEfficiencyEpochs[eyeId][telId].size() - 2;
  int lastBin = telEpochIt->second.size() - 2;

  std::vector<double> p;
  if (thisBin == -1) {
    for (unsigned int i = 0; i < 3; ++i)
      //p.push_back(fPerTelescopeOpticalEfficiencyCorrections[eyeId][telId][3*lastBin + i]);
      p.push_back(telIt->second[3*lastBin + i]);
  } else {
    for (unsigned int i = 0; i < 3; ++i)
      //p.push_back(fPerTelescopeOpticalEfficiencyCorrections[eyeId][telId][3*thisBin + i]);
      p.push_back(telIt->second[3*thisBin + i]);
  }

  /*
  cout <<"thisBin: " << thisBin <<"\n";
  cout << "gpsSecond " << gpsSecond << "\n";
  cout << p[0] << " " << p[1] << " " << p[2] <<"\n";
  */
  //unsigned int lastEpoch = fPerTelescopeOpticalEfficiencyEpochs[eyeId][telId][lastBin+1];
  unsigned int lastEpoch = telEpochIt->second[lastBin+1];
  int x;

  if (thisBin > -1)    // means that  gpsSecond is within Epoch range
    x = gpsSecond - 1016048411;
  else                 // means that  gpsSecond is greater than Epoch range
    x = lastEpoch - 1016048411;

  const double cali =
    p[0] + p[1] * x + p[2] * pow(x, 2);

  value = cali;

  return true;
}


VManager::Status
FParametricOpticalEfficiencyLossManager::GetData(vector<int>& returnData,
                                                 const string& componentProperty,
                                                 const string& componentName,
                                                 const IndexMap& /*componentIndex*/)
  const
{
  if (componentProperty != "status" || componentName != "fd_optical_efficiency")
    return eNotFound;
  returnData.clear();
  returnData.resize(441, 0);  // eGood
  return eFound;
}


VManager::Status
FParametricOpticalEfficiencyLossManager::GetData(int& returnData,
                                                 const string& componentProperty,
                                                 const string& componentName,
                                                 const IndexMap&)
  const
{
  // this manager only knows about fd_optical_efficiency
  if (componentName != "fd_optical_efficiency")
    return eNotFound;

  // check property availability
  if (fAvailableComponents.find(componentProperty) == fAvailableComponents.end())
    return eNotFound;

  // parametric optical efficiency loss correction is valid only for the given time
  const TimeStamp detTime = Detector::GetInstance().GetTime();
  const unsigned int gpsSecond = detTime.GetGPSSecond();

  if (componentProperty == "start_time") {
    returnData = gpsSecond-1;
    return eFound;
  }

  if (componentProperty == "end_time") {
    // +1 as a margin
    returnData = gpsSecond+1;
    return eFound;
  }

  return eNotFound;
}