FdSimEventChecker.cc

Go to the documentation of this file.

#include "FdSimEventChecker.h"

#include <utl/ErrorLogger.h>
#include <utl/RandomEngine.h>

#include <det/Detector.h>

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

#include <atm/InclinedAtmosphericProfile.h>
#include <atm/ProfileResult.h>

#include <utl/TimeStamp.h>
#include <utl/Reader.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.h>

#include <fevt/FEvent.h>
#include <fevt/Header.h>
#include <fevt/Eye.h>
#include <fevt/Telescope.h>
#include <fevt/EyeHeader.h>

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <fwk/RandomEngineRegistry.h>
#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <utl/UTMPoint.h>

#include <CLHEP/Random/Randomize.h>

#include <iostream>

#include <TH2.h>
#include <TFile.h>

// #define DEBUGTHIS

#ifdef DEBUGTHIS
#  include <TH2D.h>
#endif

using CLHEP::RandFlat;

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


VModule::ResultFlag
FdSimEventChecker::Init()
{
  // get data cards

  CentralConfig* const cc = CentralConfig::GetInstance();
  Branch topB = cc->GetTopBranch("FdSimEventChecker");

  if (!topB) {
    ERROR("Could not find branch FdSimEventChecker");
    return eFailure;
  }

  fVerbosity = 0;
  if (topB.GetChild("verbosity"))
    topB.GetChild("verbosity").GetData(fVerbosity);

  // laser deselection flags
  topB.GetChild("upTime").GetData(fUpTimeCheck);
  topB.GetChild("geometry").GetData(fGeometryCheck);
  topB.GetChild("geometryCherenkovHEAT").GetData(fGeometryCheckCherenkoHEAT);
  topB.GetChild("geometryCherenkovHECO").GetData(fGeometryCheckCherenkoHECO);
  topB.GetChild("fastSkip").GetData(fFastSkip);
  topB.GetChild("fastRepeat").GetData(fFastRepeat);
  topB.GetChild("useCDAS").GetData(fUseCDAS);

  hMaxVA = NULL;
  if (fGeometryCheckCherenkoHECO && topB.GetChild("maxVAfileHECO")) {
    string f;
    topB.GetChild("maxVAfileHECO").GetData(f);
    TFile* f_eff = new TFile(f.c_str());
    hMaxVA = (TH2D*)f_eff->Get("hMaxVA");
    hMaxVA->SetDirectory(0);
    f_eff->Close();
  }

  fMinViewingAngle = 1e10; // large number -> no cut
  if (topB.GetChild("minViewingAngle"))
    topB.GetChild("minViewingAngle").GetData(fMinViewingAngle);

  fRandomEngine = &RandomEngineRegistry::GetInstance().Get(RandomEngineRegistry::eDetector);


  // -------- do some output ------------
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n";
  info << "            -- shower mode -- \n"
          "              up-time check: " << fUpTimeCheck << "\n"
          "    geometry check (std FD): " << fGeometryCheck << "\n"
          "geometry check (PCGF, HEAT): " << fGeometryCheckCherenkoHEAT << "\n"
          "geometry check (PCGF, HECO): " << fGeometryCheckCherenkoHECO << "\n"
          "                   fastSkip: " << fFastSkip << " (to issue eBreakLoop) \n"
          "                 fastRepeat: " << fFastRepeat << " (to issue eContinueLoop) \n"
          "                    useCDAS: " << fUseCDAS << "\n";
  if (fMinViewingAngle<180*degree) {
    info << "    viewing angle cut: " << fMinViewingAngle/degree << " deg\n";
  }
  INFO(info);

  if (fFastSkip && fFastRepeat) {
    ERROR("It is not possible to isse at the same time eBreakLoop and eContinueLoop. Fix xml configuration.");
    return eFailure;
  }

#ifdef DEBUGTHIS
  int nE = 1000;
  int nD = 1000;
  TH2D* hDebug = new TH2D("debug", "threshold;lgE;dist/m", nE, 14.5, 17.5, nD, 0, 12000);
  for (int i=0; i<nE; ++i) {
    for (int j=0; j<nD; ++j) {
      double lgE = 14.5 + (17.5-14.5) / (nE-1) * i;
      double dist = 0. + 12000. / (nD-1) * j;
      hDebug->Fill(lgE, dist, InterpolatePCGFAcceptanceMap(lgE, dist));
    }
  }
  hDebug->SaveAs("debug.root");
#endif


  return eSuccess;
}


VModule::ResultFlag
FdSimEventChecker::Run(Event& event)
{
  if (!event.HasFEvent())
    return eContinueLoop;

  if (!event.HasSimShower()) {
    ERROR ("Event has no simulated shower.");
    return eFailure;
  }

  evt::ShowerSimData& simShower = event.GetSimShower();
  const CoordinateSystemPtr showerCS = simShower.GetShowerCoordinateSystem();
  Point showerPos(0, 0, 0, showerCS);
  Vector showerAxis(0, 0, -1, showerCS);

  bool simAny = false;
  bool skipAny = false;
  bool minViewingAngleRejected = false;

  // check if we run laser simulation
  //   ->  no geometry cutting !!!
  bool isLaser = simShower.HasLaserData();
  if (isLaser) {
    WARNING("Running in Laser simulation mode. Skipping geometry cuts.");
  }

  Detector& detector = Detector::GetInstance();
  const FDetector& detFD = detector.GetFDetector();

  if (fVerbosity >= 1) {
    DumpDB();
  }

  // --------------------------------------------------------------
  // up time check
  if (fUpTimeCheck && fUseCDAS) {

    double cdasUpTimeFraction = detFD.GetCDASUpTimeFraction();
    int cdasDaqStatus = detFD.GetCDASDAQStatus();
    double r = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 1.0);

    if (cdasDaqStatus == 0)
      cdasUpTimeFraction = 0;

    if (r > cdasUpTimeFraction) {
      ostringstream info;
      info << "Skipping event beacuse of CDAS up time cdasTimeFraction=" << cdasUpTimeFraction;
      INFO(info);
      skipAny = true;
    }
  }

  if (!skipAny) { // check for CDAS veto

    for (FDetector::EyeIterator iEye = detFD.EyesBegin() ;
         iEye != detFD.EyesEnd() ; ++iEye){

      const unsigned int eyeId = iEye->GetId();

      // --------------------------------------------------------------
      // up time check
      if (fUpTimeCheck) {
        double eyeUpTimeFraction = iEye->GetUpTimeFraction();
        int eyeDaqStatus = iEye->GetDAQStatus();
        double fdasVetoFraction = iEye->GetFDASVetoFraction();
        double cdasVetoFraction = fUseCDAS ? iEye->GetCDASVetoFraction() : 1.;
        double r = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 1.0);

        if (eyeDaqStatus == 0)
          eyeUpTimeFraction = 0;

        if (r > eyeUpTimeFraction*fdasVetoFraction*cdasVetoFraction) {
          ostringstream info;
          info << "Skipping eye=" << eyeId << " because of upTimeFraction *fdasVetoFraction * cdasVetoFraction =" << eyeUpTimeFraction*fdasVetoFraction*cdasVetoFraction;
          INFO(info);
          skipAny = true;
          continue;
        }
      }

      if (!event.GetFEvent().HasEye(eyeId, fevt::ComponentSelector::eExists)) {
        ERROR("EYE missing!");
      }
      fevt::Eye& eyeEvent = event.GetFEvent().GetEye(eyeId, fevt::ComponentSelector::eExists);
      eyeEvent.SetStatus(fevt::ComponentSelector::eInDAQ);

      // --------------------------------------------------------------
      // geometry check
      bool geometryRejected = false;
      if (fGeometryCheck && !isLaser) {

        const Point eyePosition = iEye->GetPosition();
        const double dist = (showerPos - eyePosition).GetMag();

        const double lgE = log10(simShower.GetEnergy()/eV);
        const double Rcut = Rcutoff(lgE);

        if (dist > Rcut) {

          ostringstream info;
          info << "Skipping eye=" << eyeId
               << " because of geometry cut (distance=" << dist/km << "km > R_cut=" << Rcut/km << "km)";
          INFO(info);

          geometryRejected = true;
          skipAny = true;
          eyeEvent.SetStatus(fevt::ComponentSelector::eDeSelected);
        }
      }

      if (fGeometryCheckCherenkoHEAT && !isLaser) {

        const Point eyePosition = iEye->GetPosition();
        const double dist = (showerPos - eyePosition).GetMag();

        const double lgE = log10(simShower.GetEnergy()/eV);
        const double XmaxAngle = CalculateXmaxViewingAngle(*iEye, simShower);
        const bool cut = PCGFTriggerProbabilityCut(lgE, dist, XmaxAngle);

        if (cut) {
          ostringstream info;
          info << "Skipping eye=" << eyeId
               << " because of PCGF HEAT geometry cut (distance=" << dist/km << "km, lgE=" << lgE << ", angle=" << XmaxAngle/deg << "deg)";
          INFO(info);

          geometryRejected = true;
          skipAny = true;
          eyeEvent.SetStatus(fevt::ComponentSelector::eDeSelected);
        }
      }

      if (fGeometryCheckCherenkoHECO && !isLaser) {

        const Point HEPosition = detFD.GetEye(5).GetPosition();
        const Point COPosition = detFD.GetEye(4).GetPosition();
        const double distHE = (showerPos - HEPosition).GetMag();
        const double distCO = (showerPos - COPosition).GetMag();

        const double lgE = log10(simShower.GetEnergy()/eV);
        const double XmaxAngleHE = CalculateXmaxViewingAngle(detFD.GetEye(5), simShower);
        const double XmaxAngleCO = CalculateXmaxViewingAngle(detFD.GetEye(4), simShower);
        const bool cutHE = PCGFTriggerProbabilityCutHECO(lgE, distHE, XmaxAngleHE);
        const bool cutCO = PCGFTriggerProbabilityCutHECO(lgE, distCO, XmaxAngleCO);
        const bool cut = cutHE || cutCO;

        if (cut) {
          ostringstream info;
          info << "Skipping eye=" << eyeId
               << " because of PCGF HECO geometry cut (distHE=" << distHE/km << "km, distCO=" << distCO/km << "km, lgE=" << lgE
               << ", angleHE=" << XmaxAngleHE/deg << "deg, angleCO=" << XmaxAngleCO/deg << "deg)";
          INFO(info);

          geometryRejected = true;
          skipAny = true;
          eyeEvent.SetStatus(fevt::ComponentSelector::eDeSelected);
        }
      }

      // Only telescopes with status eInDAQ in the FEvent HERE will be considered by any Modules
      // further donwstream in the sequence !
      for (fdet::Eye::TelescopeIterator iTel = iEye->TelescopesBegin();
           iTel != iEye->TelescopesEnd(); ++iTel) {

        const unsigned int telId = iTel->GetId();

        // --------------------------------------------------------------
        // up time check
        if (fUpTimeCheck) {

          double telUpTimeFraction = iTel->GetUpTimeFraction();
          int telDaqStatus = iTel->GetDAQStatus();
          double r = RandFlat::shoot(&fRandomEngine->GetEngine(), 0.0, 1.0);

          if (telDaqStatus == 0)
            telUpTimeFraction = 0;

          if (r > telUpTimeFraction) {
            ostringstream info;
            info << "Skipping eye=" << eyeId << " tel=" << telId
                 << " because of upTimeFraction=" << telUpTimeFraction;
            INFO(info);
            skipAny = true;
            continue;
          }
        }
        // --------------------------------------------------------------
        // Min Angle check
        if (fMinViewingAngle < 180*deg) {
          if (CalculateMinViewingAngle(*iTel, simShower)>fMinViewingAngle) {
            minViewingAngleRejected = true;
            skipAny = true;
            continue;
          }
        }

        if (!eyeEvent.HasTelescope(telId, fevt::ComponentSelector::eExists)) {
          ERROR("Telescope does not exists");
        }
        fevt::Telescope& tel = eyeEvent.GetTelescope(telId, fevt::ComponentSelector::eExists);
        tel.SetStatus(fevt::ComponentSelector::eInDAQ);

        ostringstream pointing;
        pointing << "sim: " << iTel->GetTelescopePointingId();
        tel.SetRawTelPointing(pointing.str());

        if (!geometryRejected) {
          simAny = true;
        }
      } // loop tel
    } // loop eyes
  } // if not CDAS veto

  if (minViewingAngleRejected) {
    ++RunController::GetInstance().GetRunData().GetNamedCounters()["FdSimEventChecker/rejectedDueToMinViewingAngle"];
    INFO("At least one telescope was rejected because of minimum viewing angle cut.");
  }

  if (!simAny) {
    INFO("No telescopes/eyes left for simulation.");
  }

  if (skipAny) {
    ++RunController::GetInstance().GetRunData().GetNamedCounters()["FdSimEventChecker/eventsWithSkippedTelescopes"];
    INFO("Some telescopes/eyes have been removed from DAQ.");
  }

  if (!simAny) {
    ++RunController::GetInstance().GetRunData().GetNamedCounters()["FdSimEventChecker/NoFdSim"];
    if (fFastSkip)
      return eBreakLoop;  // Break the sim "loop=1". Don't waste time.
    if (fFastRepeat)
      return eContinueLoop;  // sample next geometry in the sim loop
  }

  /*  if (fMinViewingAngle<180*deg) {
    const CoordinateSystemPtr referenceCS = det::Detector::GetInstance().GetReferenceCoordinateSystem();
    const CoordinateSystemPtr localCS = simShower.GetLocalCoordinateSystem();
    utl::Point  core = simShower.GetPosition();
    const double E = simShower.GetEnergy()/giga;//  convert to Gev
    const double Zenith = (-simShower.GetDirection()).GetTheta(localCS)/degree; // converted to degree
    const double Azimuth = (-simShower.GetDirection()).GetPhi(localCS)/degree; // degree
    const double core_X= core.GetX(referenceCS)/centi ; // converted to cm
    const double core_Y= core.GetY(referenceCS)/centi ; // cm

    ofstream outfile;
    outfile.open ("FedToCORSIKA.txt", std::ios_base::app);
    outfile << E << " " << Zenith <<" "<< Azimuth << " "<< core_X << " " << core_Y << " "<< endl;
    INFO(" Sim E(GeV) , zenith(deg) ,azimuth(deg) and shower core(cm) wrote to text file to fed to CORSIKA ");
    outfile.close();
  }
  */

  return eSuccess;
}


VModule::ResultFlag
FdSimEventChecker::Finish()
{
  return eSuccess;
}


void
FdSimEventChecker::DumpDB()
  const
{
  cout << " ************** dumpDB ************* " << endl;

  Detector& detector = Detector::GetInstance();
  const FDetector& detFD = detector.GetFDetector();

  cout << " ******* time = " << detector.GetTime() << " ********* " << endl;
  cout << " ******* GPS = " << detector.GetTime().GetGPSSecond() << " ********* " << endl;

  double cdasUpTimeFraction = detFD.GetCDASUpTimeFraction();
  int cdasDaqStatus = detFD.GetCDASDAQStatus();

  cout << " CDAS uptime : " << cdasUpTimeFraction
       << " status: " << cdasDaqStatus << endl;

  for (FDetector::EyeIterator iEye = detFD.EyesBegin() ;
       iEye != detFD.EyesEnd() ; ++iEye){

    unsigned int eyeId = iEye->GetId();

    double eyeUpTimeFraction = iEye->GetUpTimeFraction();
    int eyeDaqStatus = iEye->GetDAQStatus();

    cout << " eye id=" << eyeId
         << " uptime: " << eyeUpTimeFraction
         << " status: " << eyeDaqStatus << endl;

    for (fdet::Eye::TelescopeIterator iTel = iEye->TelescopesBegin();
         iTel != iEye->TelescopesEnd(); ++iTel) {

      unsigned int telId = iTel->GetId();

      double telUpTimeFraction = iTel->GetUpTimeFraction();
      int telDaqStatus = iTel->GetDAQStatus();

      cout << " tel id=" << telId
           << " uptime: " << telUpTimeFraction
           << " status: " << telDaqStatus << endl;

    } // loop tel

  } // loop eyes

}


double
FdSimEventChecker::Rcutoff(const double lgE)
  const
{
  const double x = fmax(16.5, fmin(21., lgE));

  const double p1 =  4.86267e+05;
  const double p2 = -6.72442e+04;
  const double p3 =  2.31169e+03;

  return p1 + p2*x + p3*x*x;
}

//similar to what is done for HEAT only PCGF - TH2D used explicitely, HE+CO was used in simulations to make this
//TLT trigger probability = 1e-4 defines the max allowed VA
bool
FdSimEventChecker::PCGFTriggerProbabilityCutHECO(const double lgE, const double dist, const double angle)
  const
{
  if (hMaxVA == NULL)
    return angle > 20*degree;
  int binLogE = 0;
  int binDist = 0;
  int dz = 0;
  hMaxVA->GetBinXYZ(hMaxVA->FindBin(lgE, dist), binLogE, binDist, dz);
  if (binLogE < 1 || binLogE > hMaxVA->GetNbinsX() || binDist < 1 || binDist > hMaxVA->GetNbinsY())
    return false;
  const double thresholdAngle = hMaxVA->GetBinContent(binLogE, binDist) * degree;
  return angle > thresholdAngle; // true: cut, false : no cut
}


bool
FdSimEventChecker::PCGFTriggerProbabilityCut(const double lgE, const double dist, const double angle)
  const
{
  const double thresholdAngle = InterpolatePCGFAcceptanceMap(lgE, dist) * degree;
  return angle > thresholdAngle; // true: cut, false : no cut
}

/*
  The acceptance map is a 2D histogram binned in lgE and distance/m
  around HEAT, specifying the "Xmax viewing angle" (angle between line
  of sight from HEAT to Xmax wrt. to shower axis) above which the
  (Cherenkov) trigger threshold drops below 0.001.
 */
double
FdSimEventChecker::InterpolatePCGFAcceptanceMap(const double lgE, const double dist)
  const
{
  // this cut is based on 6M CONEX simulations. Please excuse the ASCII dump of a TH2D...
  // the cut is only effective inside the margins of the study: lgE=15...17 and dist=0...10km
  const int nEbins = 5;
  const double lgEmin = 15;
  const double lgEmax = 17;
  const int nDbins = 12;
  const double dmin = 0;
  const double dmax = 10000;

  // bin content in order: iebin + idbin * Nebins

  const double angleMax[60] = {
     17.5674
    ,63.6736
    ,47.1347
    ,100
    ,100
    ,21.0435
    ,54.6757
    ,84.3739
    ,100
    ,100
    ,0
    ,33.0216
    ,66.2782
    ,100
    ,100
    ,0
    ,21.6176
    ,42.0684
    ,80.9825
    ,43.8473
    ,0
    ,0
    ,32.9672
    ,61.9826
    ,33.5471
    ,0
    ,0
    ,0
    ,40.9694
    ,6.69125
    ,0
    ,0
    ,0
    ,0
    ,5.98639
    ,0
    ,0
    ,0
    ,0
    ,49.8777
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
    ,0
  };

  const double ebin = (lgE - lgEmin) / (lgEmax - lgEmin) * nEbins;
  const double dbin = (dist/m - dmin) / (dmax - dmin) * nDbins;

  const int iebin = int(ebin); // x.0 ... x.999  ->  x
  const int idbin = int(dbin); // x.0 ... x.999  ->  x

  if (iebin >= nEbins)
    return 360; // no cut, high-e

  if (iebin < 0 || idbin < 0 || idbin >= nDbins)
    return 0; // full cut: low-e, far away

  // linear interpolate in 2D, first in lgE, then in D

  const double dE = ebin - iebin; // 0...0.999
  const double dD = dbin - idbin; // 0...0.999

  double vecTemp[nDbins]; // this is the vector in d-direction interpolated between two lgE bins
  for (int idbin = 0; idbin < nDbins; ++idbin)
    vecTemp[idbin] = angleMax[iebin + idbin * nEbins];

  if (dE < 0.5) { // dE = 0...0.5
    if (iebin > 0) {
      // towards lower lgE
      for (int idbin = 0; idbin < nDbins; ++idbin)
        vecTemp[idbin] += (dE-0.5) * (angleMax[iebin + idbin * nEbins] - angleMax[iebin-1 + idbin * nEbins]);
    }
  } else { // dE = 0.5..0.999
    // towards higher lgE
    if (iebin < nEbins-1) {
      for (int idbin = 0; idbin < nDbins; ++idbin)
        vecTemp[idbin] += (dE-0.5) * (angleMax[iebin+1 + idbin * nEbins] - angleMax[iebin + idbin * nEbins]);
    }
  }

  double result = vecTemp[idbin];

  if (dD < 0.5) { // dD = 0...0.5
    if (idbin > 0) {
      // towards lower D
      result += (dD-0.5) * (vecTemp[idbin] - vecTemp[idbin-1]);
    }
  } else { // dD = 0.5...0.99
    if (idbin < nDbins-1) {
      // towards higher D
      result += (dD-0.5) * (vecTemp[idbin+1] - vecTemp[idbin]);
    }
  }
  return result;
}


double
FdSimEventChecker::CalculateMinViewingAngle(const fdet::Telescope& tel,
                                            const evt::ShowerSimData& shower)
  const
{
  utl::Vector telescope_axis = tel.GetAxis();
  utl::Vector shower_axis = shower.GetDirection();
  utl::Point  shower_core = shower.GetPosition();
  utl::Point Eye_pos = tel.GetPosition();

  // CHECK LightAtDiaphragmSimulator FieldOfViewChecker to be smarter !!!

  double theMinAngle = 1e20;
  const int nBins = 20000;
  for (int i = 0; i < nBins; ++i) {

    double length = 20.0*meter * (i - nBins/2);

    utl::Point pSh = shower_core + length * shower_axis;

    utl::Vector view = pSh - Eye_pos;

    double theAngleFOV = acos(view*telescope_axis/view.GetMag());

    if (theAngleFOV>tel.GetCamera().GetFieldOfView()) // 23.07 degree by default
      continue;

    theMinAngle = min(theMinAngle, (acos(shower_axis*view/view.GetMag())));
  } // along shower axis

  return theMinAngle;
}


double
FdSimEventChecker::CalculateXmaxViewingAngle(const fdet::Eye& eye,
                                             const evt::ShowerSimData& shower) const
{
  utl::Point eyePos = eye.GetPosition();
  CoordinateSystemPtr localCS = fwk::LocalCoordinateSystem::Create(eyePos);

  const Point& showerCore =  shower.GetPosition();
  // const TimeStamp& timeAtCore = shower.GetTimeStamp();
  const Vector& showerAxisReal = shower.GetDirection(); // REAL SHOWER MOVEMENT
  Vector showerAxis = showerAxisReal;
  bool upward = true;
  if (showerAxis.GetTheta(localCS)>90.*deg) {
    showerAxis *= -1;  // reco-def (ALWAYS POINTING UPWARDS)
    upward = false;
  }

  const Vector eyeToCore = showerCore-eyePos;
  Vector sdp = cross(showerAxis, eyeToCore);
  // const double rp = sdp.GetMag() * (upward ? -1 : 1);
  sdp.Normalize();

  const Vector vertical(0, 0, 1, localCS);
  utl::Vector horizontalInSDP = cross(sdp, vertical);
  horizontalInSDP.Normalize();

  const double chi0 = Angle(horizontalInSDP, showerAxis) - (upward ? 180.*deg : 0.*deg);

  if (shower.HasGHParameters()) {

    const double XmaxSim = shower.GetGHParameters().GetXMax();

    double XmaxChi = 0;
    //double
    //otoa::fd::ConvertXToChi(const double x, const fdet::Telescope& detTel,
    //const Point& core, const Vector& axis, const Vector& sdp)
    {
      utl::ReferenceEllipsoid wgs84(utl::ReferenceEllipsoid::Get(utl::ReferenceEllipsoid::eWGS84));

      try {
        const double deltaX = 10.*g/cm2;
        atm::InclinedAtmosphericProfile incModel(showerCore, showerAxis, deltaX);

        // use the inclined atmosphere model
        const atm::ProfileResult& distVsSlantDepth = incModel.EvaluateDistanceVsSlantDepth();
        //const ProfileResult &heightVsSlantDepth = incModel.EvaluateHeightVsSlantDepth();

        double distance = 0;
        if (distVsSlantDepth.MinX() < XmaxSim && XmaxSim < distVsSlantDepth.MaxX())
          distance = -distVsSlantDepth.Y(XmaxSim); // DEF: p = core + distance * axis

        // const utl::Point& telPos = detTel.GetPosition();
        const utl::Point& xPoint = showerCore + distance*showerAxis;

        //CoordinateSystemPtr localCS = fwk::LocalCoordinateSystem::Create(telPos);
        // Note that the following is wrong for telescopes:
        //const CoordinateSystemPtr localCS = detTel.GetTelescopeeCoordinateSystem();
        const Vector vertical(0, 0, 1, localCS);
        const Vector horizontal = cross(sdp, vertical);

        XmaxChi = Angle((xPoint-eyePos), horizontal);

        //return angle;
      } catch (atm::InclinedAtmosphericProfile::InclinedAtmosphereModelException& e) { 
      } catch (OutOfBoundException& e) {
      }

      //return 0;
    }

    return chi0 - XmaxChi;
  }

  return -1;
}