/RayTracer.cc

Go to the documentation of this file.

#include "RayTracer.h"
#include "Filter.h"
#include "Lens.h"
#include "Mirror.h"
#include "Camera.h"

#include <fdet/Telescope.h>

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


#include <fwk/RunController.h>
#include <evt/Event.h>
#include <evt/Header.h>

#include <TPolyLine3D.h>
#include <TCanvas.h>
#include <TObjArray.h>
#include <TFile.h>
#include <TRandom.h>

#include <string>
#include <iostream>

using namespace std;
using namespace utl;
using namespace TelescopeSimulatorKG2;


void
RayTracer::Track::AddPoint(const utl::Point& p)
{
  fX.push_back(p.GetX(fCS));
  fY.push_back(p.GetY(fCS));
  fZ.push_back(p.GetZ(fCS));
}

TPolyLine3D*
RayTracer::Track::GetLine()
{
  TPolyLine3D* track = new TPolyLine3D(fX.size(), &fX.front(), &fY.front(), &fZ.front());
  track->SetLineColor(fColor);
  return track;
}


RayTracer::RayTracer(const int Verbosity,
                     const fdet::Telescope& tel,
                     utl::RandomEngine& random,
                     bool simulateShadow,
                     bool simualteCameraSupportShadow,
                     bool simulateMercedesStars,
                     bool simulateFilterStructure,
                     bool simulateHaloEffects,
                     bool simulateGhostEffects,
                     const int maxMirrorReflections,
                     const double mirrorSize,
                     const double mirrorSegmentSigma,
                     const double mirrorRadiusSigma,
                     const double mirrorAbsorptionTop,
                     const double mirrorAbsorptionBot,
                     utl::TabulatedFunction* mirrorDiffusionTop,
                     utl::TabulatedFunction* mirrorDiffusionBot,
                     const double filterIncreaseReflInside,
                     const double filterIncreaseReflOutside,
                     const double filterPosition,
                     const double filterPositionVertical,
                     const double filterPositionVertical2,
                     const double filterPositionHorizontal,
                     const double filterDustAbsorption,
                     const double lensIncreaseRefl,
                     const double lensPosition,
                     const double minLensThickness,
                     const double torusRadius,
                     const double tubeRadius,
                     const double torusZ0,
                     const double pmt_n,
                     bool plotPhotonTracks,
                     bool plotPhotonTracksAtMercedes,
                     double drawPhotonsProbabilty) :
  fVerbosity(Verbosity),
  fTel(&tel),
  fRandom(&random),
  fSimulateCameraShadow(simulateShadow),
  fSimualteCameraSupportShadow(simualteCameraSupportShadow),
  fSimulateMercedesStars(simulateMercedesStars),
  fSimulateFilterStructure(simulateFilterStructure),
  fSimulateHaloEffects(simulateHaloEffects),
  fSimulateGhostEffects(simulateGhostEffects),
  fMaxMirrorReflections(maxMirrorReflections),
  fPlotPhotonTracks(plotPhotonTracks),
  fDrawPhotonProbability(drawPhotonsProbabilty),
  fNTooManyReflections(0)
{
  fHasCorrectorRing = fTel->HasCorrectorRing();

  fFilter = new Filter(random,
                       tel,
                       filterIncreaseReflInside,
                       filterIncreaseReflOutside,
                       filterPosition,
                       filterPositionVertical,
                       filterPositionVertical2,
                       filterPositionHorizontal,
                       filterDustAbsorption,
                       fSimulateFilterStructure);

  fLens = new Lens(random,
                   tel,
                   fSimulateHaloEffects,
                   lensIncreaseRefl,
                   lensPosition,
                   minLensThickness,
                   torusRadius,
                   tubeRadius,
                   torusZ0);

  fMirror = new Mirror(random,
                       tel,
                       mirrorSize,
                       mirrorSegmentSigma,
                       mirrorRadiusSigma, 
                       mirrorAbsorptionTop,
                       mirrorAbsorptionBot,
                       mirrorDiffusionTop,
                       mirrorDiffusionBot);

  fCamera = new Camera(random,
                       tel,
                       plotPhotonTracksAtMercedes,
                       fSimulateMercedesStars,
                       fSimulateHaloEffects,
                       fSimulateGhostEffects,
                       pmt_n);

  if (plotPhotonTracks)
    fObjectsPhotons = new TObjArray(100);
  else
    fObjectsPhotons = nullptr;
}


RayTracer::~RayTracer()
{
  if (fPlotPhotonTracks) {
    ostringstream fnameBase, fnameRoot, fnameCanvas;
    const string id = fwk::RunController::GetInstance().GetCurrentEvent().GetHeader().GetId();
    fnameBase << "raytrace_" << id
              << "_Eye" << fTel->GetEyeId()
              << "_Tel" << fTel->GetId();
    fnameRoot << fnameBase.str() << ".root";
    fnameCanvas << fnameBase.str() << ".eps";
    TCanvas c("view3d");
    fCamera->Draw();
    fLens->Draw();
    fMirror->Draw();
    fFilter->Draw();
    for (int iObj = 0; iObj < fObjectsPhotons->GetEntriesFast(); ++iObj) {
      TPolyLine3D* photon = dynamic_cast<TPolyLine3D*>((*fObjectsPhotons)[iObj]);
      if (photon)
        photon->Draw();
    }
    c.Print(fnameCanvas.str().c_str());
    TFile f(fnameRoot.str().c_str(), "RECREATE");
    c.Write();
    f.Close();
  }

  delete fFilter;
  delete fLens;
  delete fMirror;
  delete fCamera;

  if (fVerbosity > 0){
    if (fNTooManyReflections) {
      ostringstream bugMsg;
      bugMsg << fNTooManyReflections
             << " photon(s) lost because of too many mirror reflections (>" << fMaxMirrorReflections << ")";
      INFO(bugMsg);
    }
  }

}

RTResult
RayTracer::Trace(const utl::Photon& photonIn,
                 utl::Photon& photon, // OUT
                 int& nMercReflections,
                 int& nBackScattered,
                 int& col, int& row)
{
  // count mirror reflections to not get trapped in infinite reflections while simulate Halo
  nMercReflections = 0;
  nBackScattered = 0; // -> halo, ghost
  //fmaxMirrorReflections (= 5 , 3 = up to "ghost")
  int nMirrorReflections = 0;

  bool plotTrack = fPlotPhotonTracks;
  int trackColor = kBlue;
  //bool plotThisTrack = false;
  if (fPlotPhotonTracks) {
    //plotThisTrack = (gRandom->Rndm() < fDrawPhotonProbability);
    const utl::CoordinateSystemPtr& telCS = fTel->GetTelescopeCoordinateSystem();
    if (photonIn.GetPosition().GetR(telCS) < 0.85 * utl::meter) {
      trackColor = kRed;
    }
  }

  // the main ray-tracing loop
  photon = photonIn;
  RTNext next = eFilter;
  bool finished = false;
  bool detected = false;
  bool terminate = false;

  // for plotting
  RayTracer::Track photonTrack(fTel->GetTelescopeCoordinateSystem());
  if (plotTrack) {
    photonTrack.SetColor(trackColor);
    photonTrack.AddPoint(photon.GetPosition());
    fCamera->SetPhotonTrack(&photonTrack);
    fLens->SetPhotonTrack(&photonTrack);
  }

  if (fVerbosity > 1)
    cout << "RayTracer ------- new photon ------- " << endl;

  int nReflect = 0; // to count multi-reflections between lens+filter

  while (!finished) {

    if (fVerbosity > 1){
      const utl::CoordinateSystemPtr& telCS = fTel->GetTelescopeCoordinateSystem();
      cout << "RayTracer: w=" << setprecision(3) << setw(6) << photon.GetWeight()
           << " x=" <<  setprecision(3) << setw(6) << photon.GetPosition().GetX(telCS)
           << " y=" <<  setprecision(3) << setw(6) << photon.GetPosition().GetY(telCS)
           << " z=" <<  setprecision(3) << setw(6) << photon.GetPosition().GetZ(telCS)
           << " nx=" <<  setprecision(3) << setw(6) << photon.GetDirection().GetX(telCS)
           << " ny=" <<  setprecision(3) << setw(6) << photon.GetDirection().GetY(telCS)
         << " nz=" <<  setprecision(3) << setw(6) << photon.GetDirection().GetZ(telCS)
           << " t=" << ios::fixed << photon.GetTime().GetInterval()
           <<  " next=" << RTNextName[next]  << endl;
    }

    // only for reflection counting:
    const utl::CoordinateSystemPtr& telCS = fTel->GetTelescopeCoordinateSystem();
    const bool backwards = photon.GetDirection().GetZ(telCS) > 0; // toward aperature

    // ----- debug
    /*cout << " Photon nMirRefl=" << nMirrorReflections <<  " nRefl=" << nReflect
         << " angle=" << ((backwards ? 1 : -1)*photon.GetDirection()).GetTheta(telCS) * 180/3.141 << " ";
*/

    switch (next) {

    case eFilter:
      next = fFilter->Trace(photon, photon);
      //if (backwards && next==eLens)
      //++nReflect;
      break;

    case eLens:
      //TraceWithTorus can also be called without halo. The difference was negligible in drum simulation.
      //For now always use TraceWithTorus
      next = fLens->TraceWithTorus(photon, photon);
      /*
      if (fSimulateHaloEffects) {
        //cout << "fLens->TraceWithTorus" << endl;
        next = fLens->TraceWithTorus(photon, photon);
        //next = fLens->Trace(photon, photon);
      } else
        next = fLens->TraceWithTorus(photon, photon);
      */
      if (!backwards && next==eFilter)
        ++nReflect;
      /*if (next == eLost)
        cout << " Photon with " << nReflect << " reflections nMirrorReflections=" << nMirrorReflections << " w=" << photon.GetWeight() << endl;*/
      break;

    case eCameraShadow:
      next = fCamera->TraceShadow(photon, photon, fSimualteCameraSupportShadow);
      break;

    case eMirror:
      next = fMirror->Trace(photon, photon);
      ++nMirrorReflections;
      break;

    case eCamera:
      next = fCamera->Trace(photon, photon, nMercReflections, nBackScattered, col,row);
      break;

    case eLost:
      finished = true;
      if (plotTrack)
        photonTrack.SetColor(kGray+1);
      break;

    case eDetected:
      detected = true;
      finished = true;
      break;

    case eDrawAndExit:
      //plotThisTrack = true;
      if (fObjectsPhotons)
        fObjectsPhotons->Clear();
      finished = true;
      terminate = true;
      break;

    case SIZE_RTNext:
      //just to silence warning do nothing
      ostringstream err;
      err << "do not use SIZE_RTNext for anything but output";
      ERROR(err);
      next = eLost;
      break;

    }; // end switch


    if (plotTrack /* && next!=eLost*/) { // maybe nicer, but confusing in the details ...
      photonTrack.AddPoint(photon.GetPosition());
    }

    if (!detected && nMirrorReflections > fMaxMirrorReflections) {
      //      cout << "too many mirror reflections: " << nMirrorReflections
      //     << ", max = " << fMaxMirrorReflections << endl;
      ++fNTooManyReflections;
      /*if (plotTrack)
        delete photonTrack;*/
      return eTooManyReflections;
    }

  } // end main raytracing loop

  // RU: debug
  //if (nReflect>2) {
    //photon.SetSource(3333);
  //}

  if (/*plotThisTrack ||*/ photon.GetSource() == 3333) {// secret number !!!

    fCamera->SetPhotonTrack(0);
    fLens->SetPhotonTrack(0);
    //if (nMirrorReflections>1)
    TPolyLine3D* track = photonTrack.GetLine();
    if (track && fObjectsPhotons) {
      fObjectsPhotons->Add(track);
    }
    //else
    //delete photonTrack;
  }

  if (terminate)
    return eTerminate;

  return detected ? eOK : eAbsorbed;

} // end of TelescopeSimulator::RayTracing