TelescopeSimulatorKG2/TelescopeSimulator.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "TelescopeSimulator.h"
#include "RayTracer.h"

#include <utl/Reader.h>
#include <utl/ErrorLogger.h>
#include <utl/AugerUnits.h>
#include <utl/MathConstants.h>
#include <utl/Math.h>
#include <utl/PhysicalConstants.h>
#include <utl/PhysicalFunctions.h>
#include <utl/Point.h>
#include <utl/TabulatedFunction.h>
#include <utl/TabulatedFunctionErrors.h>
#include <utl/Trace.h>
#include <utl/MultiTabulatedFunction.h>
#include <utl/Photon.h>
#include <utl/RandomEngine.h>
#include <utl/UTMPoint.h>

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

#include <det/Detector.h>

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

#include <evt/Event.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/TelescopeSimData.h>
#include <fevt/Telescope.h>
#include <fevt/PixelSimData.h>
#include <fevt/Pixel.h>

#include <boost/tuple/tuple.hpp>

#include <TDirectory.h>
#include <TFile.h>
#include <TTree.h>

#include <CLHEP/Random/Randomize.h>

#include <iomanip>
#include <fstream>
#include <sstream>

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

using CLHEP::RandFlat;



TelescopeSimulator::TelescopeSimulator() :
  fVerbosityLevel(0)
{
}


VModule::ResultFlag
TelescopeSimulator::Init()
{
  ostringstream info;
  info << " Version: " << GetVersionInfo(VModule::eRevisionNumber) << "\n";

#if 0
  fDrumMode = false;
#endif


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

  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("TelescopeSimulatorKG2");

  info << " *** Output format: Parameter value (default value) ***\n";

  //StoreLightComponentsAtPixels default value
  fStoreLightComponentsAtPixels = false;
  topB.GetChild("StoreLightComponentsAtPixels").GetData(fStoreLightComponentsAtPixels);
  info << " * Store light components at pixels: " << (!fStoreLightComponentsAtPixels ? "no" : "yes") << " (no)\n";

  //MirrorSize
  topB.GetChild("MirrorSize").GetData(fMirrorSize);
  info << " * Mirror size: " << (fMirrorSize/deg) << " deg (33.07)" << "\n";

  //Verbosity default value
  if (topB.GetChild("VerbosityLevel"))
    topB.GetChild("VerbosityLevel").GetData(fVerbosityLevel);
  info << " * Verbosity: " << fVerbosityLevel << " (0)\n";

  //outputFileName for shadow data
  fShadowDataOutName = "shadowData.root";
  if (topB.GetChild("ShadowDataOut"))
    topB.GetChild("ShadowDataOut").GetData(fShadowDataOutName);

  //RayTracingOptions default values
  fSimulateCameraShadow = true;
  fSimualteCameraSupportShadow = true;
  fSimulateMercedesStars = true;
  fSimulateFilterStructure = true;
  fSimulateHaloEffects = true;
  fSimulateGhostEffects = true;
  fMaxMirrorReflections = 1; // >=3 simulates ghost as well
  
  if (topB.GetChild("RayTracingOptions")){
    Branch optionsB = topB.GetChild("RayTracingOptions");
    optionsB.GetChild("SimulateCameraShadow").GetData(fSimulateCameraShadow);
    optionsB.GetChild("SimualteCameraSupportShadow").GetData(fSimualteCameraSupportShadow);
    optionsB.GetChild("SimulateMercedesStars").GetData(fSimulateMercedesStars);
    optionsB.GetChild("SimulateFilterStructure").GetData(fSimulateFilterStructure);
    optionsB.GetChild("SimulateHaloEffects").GetData(fSimulateHaloEffects);
    optionsB.GetChild("SimulateGhostEffects").GetData(fSimulateGhostEffects);
    optionsB.GetChild("MaxMirrorReflections").GetData(fMaxMirrorReflections);
  }
  info << " * Raytracing options: \n"
       << "             camera body shadow: " << (fSimulateCameraShadow ? "yes" : "no") << " (yes)\n"
       << "          camera support shadow: " << (fSimualteCameraSupportShadow ? "yes" : "no") << " (yes)\n"
       << "            mercedes collectors: " << (fSimulateMercedesStars ? "yes" : "no") << " (yes)\n"
       << "               filter structure: " << (fSimulateFilterStructure ? "yes" : "no") << " (yes)\n"
       << "                   halo effects: " << (fSimulateHaloEffects ? "yes" : "no") << " (yes)\n"
       << "                  ghost effects: " << (fSimulateGhostEffects ? "yes" : "no") << " (yes)\n"
       << "   max nr of mirror reflections: " << fMaxMirrorReflections << " (1)\n";
  
  
  // ---------------------------------- ** HALO INPUT ** ----------------------------------

  //MirrorOptions default values
  //fMirrorDiffusion != NULL needed for halo
  //for other parameters: values other than default only needed for abnormal modifications of the halo
  fMirrorSegmentSigma = 0.;
  fMirrorRadiusSigma = 0.;
  fMirrorAbsorptionTop = 0.;
  fMirrorAbsorptionBot = 0.;
  fMirrorDiffusionTop = nullptr;
  fMirrorDiffusionBot = nullptr;

  if (fSimulateHaloEffects) {
    if (topB.GetChild("MirrorOptions")){
      Branch mirrorOptionsB = topB.GetChild("MirrorOptions");
      
      mirrorOptionsB.GetChild("MirrorSegmentSigma").GetData(fMirrorSegmentSigma);
      mirrorOptionsB.GetChild("MirrorRadiusSigma").GetData(fMirrorRadiusSigma);
      mirrorOptionsB.GetChild("MirrorAbsorptionTop").GetData(fMirrorAbsorptionTop);
      mirrorOptionsB.GetChild("MirrorAbsorptionBot").GetData(fMirrorAbsorptionBot);
      
      //only if MirrorDiffusionTop/Bot branch exists get values for mirror diffusion
      Branch mirrorDiffusionTopB = mirrorOptionsB.GetChild("MirrorDiffusionTop");
      mirrorDiffusionTopB.GetChild("MirrorDiffusionValues").GetData(fMirrorDiffusionTopValues);
      mirrorDiffusionTopB.GetChild("MirrorDiffusionAngles").GetData(fMirrorDiffusionTopAngles);
      fMirrorDiffusionTop = new TabulatedFunction(fMirrorDiffusionTopValues, fMirrorDiffusionTopAngles);
      
      Branch mirrorDiffusionBotB = mirrorOptionsB.GetChild("MirrorDiffusionBot");
      mirrorDiffusionBotB.GetChild("MirrorDiffusionValues").GetData(fMirrorDiffusionBotValues);
      mirrorDiffusionBotB.GetChild("MirrorDiffusionAngles").GetData(fMirrorDiffusionBotAngles);
      fMirrorDiffusionBot = new TabulatedFunction(fMirrorDiffusionBotValues, fMirrorDiffusionBotAngles);
    } else {
      ostringstream err;
      err << "fMirrorDiffusion != 0 needed for halo. \n"
          << "Can be set in the TelescopeSimulator xml file. \n"
          << "Stop Simulation.";
      ERROR(err);
      return eFailure;
    }
  } //end of if fSimulateHaloEffects
  
  //Mirror output
  info << " * Mirror options:" << "\n"
       << "       segment sigma: " << fMirrorSegmentSigma / deg << " deg" << " (0.0)" << "\n"
       << "        radius sigma: " << fMirrorRadiusSigma / mm << " mm" << " (0.0)" << "\n"
       << "      absorption top: " << 100 * fMirrorAbsorptionTop << "\% (0.0)\n"
       << "   absorption bottom: " << 100 * fMirrorAbsorptionBot << "\% (0.0)\n";

  info <<  "   mirror diffusion top function: ";
  if (fMirrorDiffusionTop == nullptr)
    info << "no" << "\n";
  else 
    info << "yes" << "\n";
      
  info <<  "   mirror diffusion bottom function: ";
  if (fMirrorDiffusionBot == nullptr)
    info << "no" << "\n";
  else 
    info << "yes" << "\n";

  //FilterOptions default values
  //Might change in decision of what is standard mode (probably without halo ???)
  fFilterPosition = 12. * cm;
  fFilterPositionVertical = 1 * cm;
  fFilterPositionVertical2 = 1 * cm;
  fFilterPositionHorizontal = 1 * cm;
  fFilterIncreaseReflectionInside = 1.;
  fFilterIncreaseReflectionOutside = 1.;
  fFilterDustAbsorptionOutside = 0.;

  if (fSimulateHaloEffects) {
    if (topB.GetChild("FilterOptions")){
      Branch filterOptionsB = topB.GetChild("FilterOptions");
      filterOptionsB.GetChild("FilterIncreaseReflectionInside").GetData(fFilterIncreaseReflectionInside);
      filterOptionsB.GetChild("FilterIncreaseReflectionOutside").GetData(fFilterIncreaseReflectionOutside);
      filterOptionsB.GetChild("FilterPosition").GetData(fFilterPosition);
      filterOptionsB.GetChild("FilterPositionVertical").GetData(fFilterPositionVertical);
      filterOptionsB.GetChild("FilterPositionVertical2").GetData(fFilterPositionVertical2);
      filterOptionsB.GetChild("FilterPositionHorizontal").GetData(fFilterPositionHorizontal);
      filterOptionsB.GetChild("FilterDustAbsorptionOutside").GetData(fFilterDustAbsorptionOutside);
    }
  } //end of if fDoNoHalo

  info << " * Filter options:" << "\n"
       << "      increased reflectivity inside: " << fFilterIncreaseReflectionInside << " (1.0)" << "\n"
       << "     increased reflectivity outside: " << fFilterIncreaseReflectionOutside << " (1.0)" << "\n"
       << "             position (z-direction): " << fFilterPosition / cm << " cm" << " (12.0)" << "\n"
       << "              position vertical bar: " << fFilterPositionVertical / cm << " cm" << "\n"
       << "             position vertical bar2: " << fFilterPositionVertical2 / cm << " cm" << "\n"
       << "            position horizontal bar: " << fFilterPositionHorizontal / cm << " cm" << "\n"
       << "   absorption on dust layer outside: " << fFilterDustAbsorptionOutside * 100 << "% (0.0)" << "\n";

  //LensOptions default values
  fLensIncreaseReflection = 1.;
  fLensPosition = 0. * mm;
  fMinLensThickness = 6.14 * mm;
  fTorusRadius = 795.27 * mm;
  fTubeRadius = 8383.00 * mm;
  fTorusZ0 = fLensPosition + fTubeRadius + 5.96 * mm;
  
  if (fSimulateHaloEffects) {
    if (topB.GetChild("LensOptions")){
      Branch lensOptionsB = topB.GetChild("LensOptions");
      lensOptionsB.GetChild("LensIncreaseReflection").GetData(fLensIncreaseReflection);
      lensOptionsB.GetChild("LensPosition").GetData(fLensPosition);

      Branch torusOptionsB = lensOptionsB.GetChild("TorusParameters");
      torusOptionsB.GetChild("MinLensThickness").GetData(fMinLensThickness);
      torusOptionsB.GetChild("TorusRadius").GetData(fTorusRadius);
      torusOptionsB.GetChild("TubeRadius").GetData(fTubeRadius);
      torusOptionsB.GetChild("TorusZ0").GetData(fTorusZ0);
    }
  }
  info << " * Lens options:" << "\n"
       << "      increase reflectivity: " << fLensIncreaseReflection << " (1.0)" << "\n"
       << "     position (z-direction): " << fLensPosition / mm << " mm" << " (0.0)" << "\n"
       << "              min thickness: " << fMinLensThickness / mm << " mm" << " (6.14)" << "\n"
       << "               torus radius: " << fTorusRadius / mm << " mm" << " (795.27)" << "\n"
       << "                tube radius: " << fTubeRadius / mm << " mm" << " (8383.00)" << "\n"
       << "   position of torus center: " << fTorusZ0 / mm << " mm" << " (8388.96)" << "\n";

  //CameraOptions default values
  fPMT_n = 2.58;
  
  if (fSimulateHaloEffects) {
    if (topB.GetChild("CameraOptions")){
      Branch cameraOptionsB = topB.GetChild("CameraOptions");
      cameraOptionsB.GetChild("pmt_n").GetData(fPMT_n);
    }
  }
  info << " * Camera options: \n"
       << "   PMT window refractive index: " << fPMT_n << " (2.58)\n";

  // ---------------------------------- ** END OF HALO INPUT ** ----------------------------------
  
  //DrawingOptions default values
  fDrawNumberOfPhotons = 0;
  fDraw = false;
  fDrawMercReflections = false;
  
  if (topB.GetChild("DrawingOptions")) {
    Branch drawingB = topB.GetChild("DrawingOptions");
    fDraw = true;
    drawingB.GetChild("mercedesReflections").GetData(fDrawMercReflections);
    drawingB.GetChild("numberOfPhotons").GetData(fDrawNumberOfPhotons);
  }

  if (fDraw) {
  info << " Drawing options: \n"
            "   plot mercedes reflections: " << (!fDrawMercReflections ? "no" : "yes") << " (no)\n"
            "           save 3D plot with: " << fDrawNumberOfPhotons << " (0) raytraced photons\n";
  }
  
  INFO(info);

  // ---------------------------------------------------
  // the model config signature
  ostringstream configSignature;
  // use configSignature of normal TelescopeSimulatorKG at the moment
  // to have its drum calibration constants available
  // TO DO: enter halo here, once calibration constants are available
#warning TO DO: all changeable values have to enter all in configSignature of TelescopeSimulator.cc

  configSignature << "TelescopeSimulatorKG2";
#if 0 
  if (fDoNoShadow || fDoNoShadowSupport || fHasNoMercedes || !fDoNoFilterStructure || !fDoNoHalo || !fDoNoGhost)
    configSignature << "(camera-shadow / support-shadow / mercedes / filter-structure / halo / ghost) = ("
                    << (!fDoNoShadow ? "yes" : "no") << " / "
                    << (!fDoNoShadowSupport ? "yes" : "no") << " / "
                    << (!fHasNoMercedes ? "yes" : "no") << " / "
                    << (!fDoNoFilterStructure ? "yes" : "no") << " / "
                    << (!fDoNoHalo ? "yes" : "no") << " / "
                    << (!fDoNoGhost ? "yes" : "no") << ")";
#endif
  fGeneralConfigSignature = configSignature.str();
  INFO(configSignature);

  return eSuccess;
}


VModule::ResultFlag
TelescopeSimulator::Run(evt::Event& event)
{
  if (!event.HasFEvent()) {
    ERROR("Event has no FEvent.");
    return eFailure;
  }
  fevt::FEvent& fEvent = event.GetFEvent();

#if 0
  if (!fDrumMode && !event.HasSimShower()) {
    /*
      In drum mode (AND SPOT MODE) the handling of photons is different from normal mode.
      While normaly photons are tracked with a weight, in drum mode photons are tracked
      as physical objects that can only get absorbed - or not.
     */
    fDrumMode = true;
    cout << "fDrumMode = true! " << endl;
  }
#endif

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

  // loop eyes
  for (fevt::FEvent::EyeIterator iEye = fEvent.EyesBegin(fevt::ComponentSelector::eInDAQ);
       iEye != fEvent.EyesEnd(fevt::ComponentSelector::eInDAQ); ++iEye) {

    if (iEye->GetStatus() == fevt::ComponentSelector::eDeSelected)
      continue;

    const unsigned int eyeId = iEye->GetId();
    fevt::Eye& eyeEvent = *iEye;
    //const fdet::Eye& detEye = detFD.GetEye(*iEye);

    // loop telescopes
    for (fevt::Eye::TelescopeIterator iTel = eyeEvent.TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         iTel != eyeEvent.TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {

      const unsigned int telId = iTel->GetId();
      fevt::Telescope& telEvent = *iTel;
      const fdet::Telescope& detTel = detFD.GetTelescope(telEvent);

      const TabulatedFunction& telMeasEfficiency  = detTel.GetMeasuredRelativeEfficiency();

      fevt::TelescopeSimData& telSim = telEvent.GetSimData();
      telSim.SetConfigSignatureStr(detTel.GetConfigSignatureStr(fGeneralConfigSignature));

      const double normWavelength = detFD.GetReferenceLambda();
      //const TabulatedFunction& mirRef = detTel.GetMirror().GetReflectivity();
      const TabulatedFunction& filtTrans = detTel.GetFilter().GetTransmittance();

      // corrector

      const fdet::Camera& detCamera = detTel.GetCamera();
      const unsigned int telTraceNBins = telSim.GetNumberOfPhotonBins();
      const double telTraceBinWidth = detCamera.GetFADCBinSize();

      const double nPhotons = telSim.GetNPhotons();

      if (!nPhotons)
        continue;

      ostringstream info;
      /*info << "Raytracing eye=" << eyeId
           << ", telescope=" << telId
           << ", rDia=" << setw(6) << detTel.GetDiaphragmRadius()
           << ", nPhotons=" << nPhotons;*/
      info << "Raytracing telescope with " << nPhotons << " Photons";
      INFO(info);

      const double minWl = detFD.GetModelMinWavelength();
      const double maxWl = detFD.GetModelMaxWavelength();

      map<int, int> rayTracerResults; // for verbose info output: [status:counts]

      const double drawPhotonProbablility = (double)fDrawNumberOfPhotons / nPhotons;
      
      RayTracer* rayTracer = nullptr;
      
      const int index = telId + eyeId * 100;
      fNHit[index] = 0;
      fN_1[index] = 0;
      fN_2[index] = 0;
      fN_3[index] = 0;
      fN_4[index] = 0;
      fN_5[index] = 0;
      fN_6[index] = 0;


      unsigned int countWlBoundary = 0;
      /*
      //These can not be accessed anymore as Trace functions of components no longer give RTResult back but RTNext remove?
      unsigned int countFilter = 0;
      unsigned int countMirror = 0;
      */
      unsigned int countPhotons = 0;

      map<int, TraceD*> traces;
      for (fevt::TelescopeSimData::PhotonIterator iPhoton = telSim.PhotonsBegin();
           iPhoton != telSim.PhotonsEnd(); ++iPhoton) {

        #warning "JD: The number of photons after which rayTracer gets refreshed needs to be optimised for performance vs geometry creeping"
        if (!(countPhotons % 100000)) {
          delete rayTracer;
          rayTracer = new RayTracer(fVerbosityLevel,
                                    detTel, *fRandomEngine,
                                    fSimulateCameraShadow, fSimualteCameraSupportShadow, fSimulateMercedesStars, 
                                    fSimulateFilterStructure, fSimulateHaloEffects, fSimulateGhostEffects,
                                    fMaxMirrorReflections,
                                    fMirrorSize, fMirrorSegmentSigma, fMirrorRadiusSigma, 
                                    fMirrorAbsorptionTop, fMirrorAbsorptionBot, fMirrorDiffusionTop, fMirrorDiffusionBot,
                                    fFilterIncreaseReflectionInside, fFilterIncreaseReflectionOutside, 
                                    fFilterPosition, 
                                    fFilterPositionVertical,
                                    fFilterPositionVertical2, 
                                    fFilterPositionHorizontal, 
                                    fFilterDustAbsorptionOutside, 
                                    fLensIncreaseReflection, fLensPosition,
                                    fMinLensThickness, fTorusRadius, fTubeRadius, fTorusZ0,
                                    fPMT_n,
                                    fDraw, fDrawMercReflections, drawPhotonProbablility);  
        }
        ++countPhotons;

        const double wavelength = iPhoton->GetWavelength();

        if (wavelength < minWl || wavelength > maxWl) {
          ++countWlBoundary;
          continue;
        }


        utl::Photon photonOut;
        int col = 0;
        int row = 0;

        int nReflections = 0; // mercedes reflections
        int nHaloReflection = 0;
        const RTResult status = rayTracer->Trace(*iPhoton, photonOut, nReflections, nHaloReflection, col, row);
        const double weight = photonOut.GetWeight();
        // DEBUG drum mode with halo
        //cout << "weight after rayTracer = " << weight << endl;

        if (status == eTerminate) {
          WARNING("********************** TERMINATED BY RAYTRACER ***************************");
          return eSuccess;
        }

        if (status != eOK ) { // photon lost
          rayTracerResults[status]++;
          continue;
        } else {
          if (status == eOK && weight == 0) {
            rayTracerResults[eAbsorbed]++;
            continue;
          }
        }

        /*
        if (nHaloReflection>0) {
          cout << "RU: nHaloReflection = " << nHaloReflection << endl;
        }
        */

        rayTracerResults[status]++;

        const unsigned int pxlId = (col - 1) * detTel.GetLastRow() + row;
        const fdet::Pixel& detPix = detTel.GetPixel(pxlId);
        const TabulatedFunction& qEff = detPix.GetQEfficiency();

        // !! this factor is to achive consistency to reconstruction !!
        // do this also for drum mode
        double efficiencyCorrection = 0;
      //  if (!fDrumMode) {
          const double modelRelEff = detTel.GetModelRelativeEfficiency(wavelength);
          if (modelRelEff)
            efficiencyCorrection = telMeasEfficiency.Y(wavelength) / modelRelEff;
      //  }
        // DEBUG drum mode with halo
        /*cout << "efficiencyCorrection = " << efffactoriciencyCorrection 
             << ", measured efficiency = " << telMeasEfficiency.Y(wavelength)
             << ", wavelength = " << wavelength;
        if (efficiencyCorrection != 0)
          cout << ", modelRelEff = " << telMeasEfficiency.Y(wavelength)/efficiencyCorrection << endl;
        else 
          cout << endl;
        */

        const double time = photonOut.GetTime().GetInterval();
        const int bin = int(time / telTraceBinWidth);

        if (fVerbosityLevel > 2) {
          cerr << " bin=" << setw(4) << bin
               << " time=" << setw(6) << time/ns
               << " weight=" << setw(4) << weight
               << " wl=" << setw(4) << wavelength/nanometer
               << " col=" << setw(3) << col
               << " row=" << setw(3) << row
               << " pixelId=" << setw(3) << pxlId
               << " telTraceBinWidth=" << setw(3) << telTraceBinWidth
               << " telTraceNBins=" << setw(5) << telTraceNBins
               << " EQ=" << qEff.Y(wavelength) / qEff.Y(normWavelength)
               << " EC=" << efficiencyCorrection
               << endl;
        }

        // check for index out of bound
        if (bin < 0 || bin >= int(telTraceNBins)) {
          if (fVerbosityLevel > 2) {
            ostringstream err;
            err << "Bin out of range. bin=" << bin
                << " weight=" << weight;
            ERROR(err);
          }
          continue;
        }

        double filterFactor = filtTrans.Y(wavelength);
        // in the new simulation, the reflections on both surfaces of the filter
        // (air-glass and glass-air) are already taken into account
        // using the refractive index of the filter material n = 1.526, the transmission
        // factor due to the surface reflections is
        // T = T_in * T_out = 0.9566 * 0.9566 = 0.915
        // Careful: n_filter = 1.526 is hard-coded in Filter.cc
        filterFactor /= 0.915;

        //PMT refractive index is by default set to 2.58, corresponding to a reflectivity R=19.5%
        //nominal refractive index given by the manufacturer is 1.52, R=4.3%
        //photons are reflected off the PMT during the ray tracing if simulation of halo is enabled (default),
        //but the loss of photons due to reflection is included in the QE measurement as well
        //to avoid double-counting reflections, we correct by pmtReflectionsCorr = 1.0/T(fPMT_n=2.58) = 1.24
        //we assume QE was measured for normal incidence, hence T is also calculated for normal incidence

        double pmtRefl = Sqr((1.0 - fPMT_n) / (1.0 + fPMT_n));
        double pmtReflectionsCorr = 1.0 / (1.0 - pmtRefl);

        double photonWeight = (weight * filterFactor * pmtReflectionsCorr * qEff.Y(wavelength) / qEff.Y(normWavelength) *
                               efficiencyCorrection);

        // DEBUG drum mode with halo
        //cout << "modified photon weight = " << photonWeight << endl;
        //cout << "filterFactor_reduced = " << filterFactor << "  ";
        //cout << "qEff: " << qEff.Y(wavelength) << " / " << qEff.Y(normWavelength)
        //     << " = " << qEff.Y(wavelength)/qEff.Y(normWavelength) << endl;

        // iPhoton is initial photon before ray tracing, photonWeight is weight after ray tracing and corrections
        // In drum mode there is no weight. Photons can only get absorbed or otherwise have weight 1.
        //

#if 0
        // TODO maybe get rid of that and make no difference anymore between drum mode and normal simulation
        //if (fDrumMode) {
          if (RandFlat::shoot(&fRandomEngine->GetEngine(), 0., 1.) > photonWeight / iPhoton->GetWeight()) {
            rayTracerResults[eAbsorbed]++;
            continue; // photon absorbed
          }
          else
            photonWeight = 1.;
        //}
#endif

        // counts total photons (using weights) at the several mercedes surfaces
        switch (nReflections) {
        case 0: fNHit[index] += photonWeight; break; // direct hit
        case 1: fN_1[index] += photonWeight; break;
        case 2: fN_2[index] += photonWeight; break;
        case 3: fN_3[index] += photonWeight; break;
        case 4: fN_4[index] += photonWeight; break;
        case 5: fN_5[index] += photonWeight; break;
        case 6: fN_6[index] += photonWeight; break;
        }

        // create pixels trace if not there yet
        if (!traces.count(pxlId)) {
          if (!telEvent.HasPixel(pxlId))
            telEvent.MakePixel(pxlId);
          fevt::Pixel& pixData = telEvent.GetPixel(pxlId);
          if (!pixData.HasSimData())
            pixData.MakeSimData();
          fevt::PixelSimData& pixSimData = pixData.GetSimData();
          if (!pixSimData.HasPhotonTrace(fevt::FdConstants::eTotal))
            pixSimData.MakePhotonTrace(telTraceNBins, telTraceBinWidth, fevt::FdConstants::eTotal);
          traces[pxlId] = &(pixSimData.GetPhotonTrace(fevt::FdConstants::eTotal));
        }
        // add to pixel's trace
        (*(traces[pxlId]))[bin] += photonWeight;

        if (fStoreLightComponentsAtPixels) {
          const fevt::FdConstants::LightSource lightSource =
            (fevt::FdConstants::LightSource)iPhoton->GetSource();
          fevt::PixelSimData& pixSimData = telEvent.GetPixel(pxlId).GetSimData();
          if (!pixSimData.HasPhotonTrace(lightSource))
            pixSimData.MakePhotonTrace(telTraceNBins, telTraceBinWidth, lightSource);
          pixSimData.GetPhotonTrace(lightSource)[bin] += photonWeight;
          if (!pixSimData.HasPhotonWeightSquareTrace(lightSource))
            pixSimData.MakePhotonWeightSquareTrace(telTraceNBins, telTraceBinWidth, lightSource);
          pixSimData.GetPhotonWeightSquareTrace(lightSource)[bin] += photonWeight*photonWeight;
        }

      } // loop Photons
      delete rayTracer;

      if (fVerbosityLevel > 0) {
        int debugTotal = 0;
        for(map<int, int>::const_iterator iRTR = rayTracerResults.begin();
            iRTR != rayTracerResults.end(); ++iRTR) {
          cout << " RayTraceResult \"" << RTResultName[iRTR->first] << "\" occured "
               << iRTR->second << " times " << endl;
          debugTotal += iRTR->second;
        }
        cout << " RayTraceResult \"total\" number of photons is " << debugTotal << endl;
        cout << " RayTraceResult wl-boundaries: " << countWlBoundary << endl;
        //cout << " RayTraceResult filter: " << countFilter << endl;
        //cout << " RayTraceResult mirror: " << countMirror << endl;
        cout << " RayTraceResult grand-total: " << debugTotal + /*countFilter + countMirror +*/ countWlBoundary << endl;
      }


      if (fVerbosityLevel > 2) {
        // Output pixels with signal
        for (map<int,TraceD*>::iterator iTrace = traces.begin();
             iTrace != traces.end(); ++iTrace) {
          const int pxlId = iTrace->first;
          const int col = (pxlId - 1) / detTel.GetLastRow() + 1; // (col-1)*lastRow + row;
          const int row = (pxlId - 1) % detTel.GetLastRow() + 1;

          cerr << "eyeId=" << eyeId << " telId=" << telId << " pixelId=" << pxlId
               << " col=" << col << " row=" << row
               << " trace:"
               << flush;

          double signal = 0;
          for (unsigned int i = 0; i < telTraceNBins; ++i) {
            if (fVerbosityLevel > 3) {
              if ((*(iTrace->second))[i])
                cerr << " i: " << i
                     << " val: " << (*(iTrace->second))[i]
                     << " |"
                     << flush;
            }
            signal += (*(iTrace->second))[i];
          }
          cerr << " total signal " << signal << endl;
        }
      } // end verbose

    } // End loop over Telescopes

  } // End loop over Eyes

  return eSuccess;
} // end of Run


VModule::ResultFlag
TelescopeSimulator::Finish()
{
  // -- shadow factor calculation --
  //if (fDrumMode) {
    //TDirectory* save = gDirectory;
    //TFile output(fShadowDataOutName.c_str(), "UPDATE");
    for (map<int, unsigned int>::const_iterator iter = fNHit.begin();
         iter != fNHit.end(); ++iter)
      {
        const int index = iter->first;
/*
        ostringstream name;
        name << "shadow_" << index;
        TTree* tree = (TTree*) output.Get(name.str().c_str());
        if (!tree) {
          tree = new TTree(name.str().c_str(), "shadow info");
          tree->Branch("direct", &fNHit[index], "direct/i");
          tree->Branch("refl1", &fN_1[index], "refl1/i");
          tree->Branch("refl2", &fN_2[index], "refl2/i");
          tree->Branch("refl3", &fN_3[index], "refl3/i");
          tree->Branch("refl4", &fN_4[index], "refl4/i");
          tree->Branch("refl5", &fN_5[index], "refl5/i");
          tree->Branch("refl6", &fN_6[index], "refl6/i");
        } else {
          tree->SetBranchAddress("direct", &fNHit[index]);
          tree->SetBranchAddress("refl1", &fN_1[index]);
          tree->SetBranchAddress("refl2", &fN_2[index]);
          tree->SetBranchAddress("refl3", &fN_3[index]);
          tree->SetBranchAddress("refl4", &fN_4[index]);
          tree->SetBranchAddress("refl5", &fN_5[index]);
          tree->SetBranchAddress("refl6", &fN_6[index]);
        }
        tree->Fill();
        tree->Write();
*/
        ostringstream info;
        info << "\n\n ------------- SHADOW INFO (telId=" << index << ") -----------" << "\n"
             << " Direct hit of focal surface: " << fNHit[index] << "\n"
             << " 1 mercedes refl.           : " << fN_1[index] << "\n"
             << " 2 mercedes refl.           : " << fN_2[index] << "\n"
             << " 3 mercedes refl.           : " << fN_3[index] << "\n"
             << " 4 mercedes refl.           : " << fN_4[index] << "\n"
             << " 5 mercedes refl.           : " << fN_5[index] << "\n"
             << " 6 mercedes refl.           : " << fN_6[index] << "\n";

        const double nTot = ((double)fN_1[index] +
                             (double)fN_2[index] +
                             (double)fN_3[index] +
                             (double)fN_4[index] +
                             (double)fN_5[index] +
                             (double)fN_6[index]);
        const double nMercedes = ((double)fN_1[index] * 1. +
                                  (double)fN_2[index] * 2. +
                                  (double)fN_3[index] * 3. +
                                  (double)fN_4[index] * 4. +
                                  (double)fN_5[index] * 5. +
                                  (double)fN_6[index] * 6.);
        const double meanRefls = (double(nMercedes) / (fNHit[index] + nTot));
        info << " Mean number of reflections from mercedes: " << meanRefls << " at tel=" << index << "\n";
        info << " -----------------------------------------------------------------";
        INFO(info);
    //  }
    //output.Close();
    //gDirectory = save;
  }
  return eSuccess;
}


#if 0
// --------------------------
// Wed Jul 26 10:51:33 CEST 2006
// this is for SG spot tables !!!!!!!!!
void
TelescopeSimulator::TransformToLocalCameraCoordinates(const double laz, const double lze, double& caz, double& cze)
{
  const double phi0 = 16.*kPi/180.;

  const double llaz = lze + 0.5*kPi;
  const double llze = laz + 0.5*kPi;

  const double xx = sin(llaz) * cos(llze);
  const double yy = sin(llaz) * sin(llze);
  const double zz = cos(llaz);

  const double x = xx;
  const double y = yy * cos(phi0) - zz * sin(phi0);
  const double z = yy * sin(phi0) + zz * cos(phi0);

  caz = asin(x);
  cze = phi0 - atan(z/y);
}
#endif