SpotPhotonGenerator.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "SpotPhotonGenerator.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 <utl/Photon.h>

#include <fwk/CentralConfig.h>
#include <fwk/RunController.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/Channel.h>
#include <fdet/Mirror.h>
#include <fdet/Filter.h>
#include <fdet/Corrector.h>

#include <evt/Event.h>
#include <evt/ShowerSimData.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 <atm/ProfileResult.h>
#include <atm/InclinedAtmosphericProfile.h>

#include <CLHEP/Random/Randomize.h>

#include <boost/tuple/tuple.hpp>

#include <TH2D.h>
#include <TDirectory.h>
#include <TFile.h>
#include <TCanvas.h>
#include <TStyle.h>

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

using namespace SpotPhotonGeneratorOG;
using namespace std;
using namespace utl;
using namespace evt;
using namespace fwk;
using namespace det;
using namespace fdet;
using namespace atm;

using CLHEP::RandFlat;


// -----------------------------------------------------------------------
VModule::ResultFlag
SpotPhotonGenerator::Init()
{
  Branch topB =
    CentralConfig::GetInstance()->GetTopBranch("SpotPhotonGenerator");

  if (!topB) {
    ostringstream err;
    err << "Missing configuration file for SpotPhotonGenerator";
    ERROR(err);
    return eFailure;
  }

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

  //
  // Initialize the spot simulation mode
  //

  topB.GetChild("NSpotPhotons").GetData(fNSpotPhotons);
  topB.GetChild("MinRDiaPhotons").GetData(fRDiaMin);
  topB.GetChild("MaxBin").GetData(fMaxBin);
  topB.GetChild("MinBin").GetData(fMinBin);

  // ------------------------------------------------
  //
  // do some final output

  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "    number of photons: " << fNSpotPhotons << "\n"
          "           min radius: " << setw(3) << fRDiaMin/m << "m\n"
          "              min bin: " << fMinBin << "\n"
          "              max bin: " << fMaxBin;
  INFO(info);

  return eSuccess;
}


// -----------------------------------------------------------------------
VModule::ResultFlag
SpotPhotonGenerator::Run(evt::Event& event)
{
  Detector& detector = Detector::GetInstance();
  const FDetector& detFD = detector.GetFDetector();

  if (!event.HasFEvent()) {  // initialize event
    ERROR("No FEvent. Check your ModuleSequence.");
    return eBreakLoop;
  }

  fevt::FEvent& fEvent = event.GetFEvent();

  static bool initialized = false;
  if (!initialized) {
    initialized = true;

    const TimeStamp& simTime = detector.GetTime();

    // set the FD event times
    fEvent.GetHeader().SetTime(simTime);
    event.GetHeader().SetTime(simTime);
    ostringstream idStr;
    //idStr << "SpotSimulation_" << fEventCounter;
    event.GetHeader().SetId(idStr.str());

    for (FDetector::EyeIterator iEye = detFD.EyesBegin(); iEye != detFD.EyesEnd(); ++iEye) {
      const unsigned int eyeId = iEye->GetId();
      fEvent.MakeEye(eyeId, fevt::ComponentSelector::eInDAQ);

      fevt::Eye& eyeEvent = fEvent.GetEye(eyeId, fevt::ComponentSelector::eInDAQ);
      for (fdet::Eye::TelescopeIterator iTel = iEye->TelescopesBegin();
           iTel != iEye->TelescopesEnd(); ++iTel) {
        const unsigned int telId = iTel->GetId();

        if (eyeEvent.HasTelescope(telId, fevt::ComponentSelector::eInDAQ))
          ERROR("Telescope already there");

        eyeEvent.MakeTelescope(telId, fevt::ComponentSelector::eInDAQ);
        fevt::Telescope& telEvent = eyeEvent.GetTelescope(telId, fevt::ComponentSelector::eInDAQ);
        telEvent.MakeSimData();

        fevt::TelescopeSimData& telSim = telEvent.GetSimData();
        telSim.SetPhotonsStartTime(simTime);
        telSim.SetNumberOfPhotonBins(fMaxBin);

      }
    }
  } // end creating FEvent

  if (DoSpot(event))
    return eSuccess;

  return eBreakLoop;
}


// -----------------------------------------------------------------------
bool
SpotPhotonGenerator::DoSpot(evt::Event& event)
{
  fevt::FEvent& fEvent = event.GetFEvent();
  Detector& detector = Detector::GetInstance();
  const FDetector& detFD = detector.GetFDetector();

  static bool finishedOneSpot = false;
  if (finishedOneSpot) {
    FinishSpot();
    // clean up
    fevt::TelescopeSimData& telSim =
      fEvent.GetEye(fCurrentEyeId,
                    fevt::ComponentSelector::eInDAQ).GetTelescope(fCurrentTelId,
                                                                  fevt::ComponentSelector::eInDAQ).GetSimData();
    telSim.ClearPhotons();
    goto innerLoop; // jump into the innermost loop to continue with spot list
  }

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

    static fevt::Eye::TelescopeIterator iTel;
    for (iTel = iEye->TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         iTel != iEye->TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {

      static unsigned int iPix;
      for ( iPix = detFD.GetTelescope(*iTel).GetFirstPixelId() ;
            iPix <= detFD.GetTelescope(*iTel).GetLastPixelId(); ++iPix) {
        { // need extra scope here

          const fdet::Telescope& detTel = detFD.GetTelescope(*iTel);
          const fdet::Pixel& detPixel = detTel.GetPixel(iPix);

          fCurrentEyeId = iEye->GetId();
          fCurrentTelId = iTel->GetId();
          fCurrentPixId = iPix;

          fevt::TelescopeSimData& telSim = iTel->GetSimData();
          TelescopeSimulatorKG::TelescopeSimulator* const raytraceModule =
            dynamic_cast<TelescopeSimulatorKG::TelescopeSimulator*>(&RunController::GetInstance().GetModule("TelescopeSimulatorKG"));
          fSpotPhotonList.clear();
          raytraceModule->SetSpotPhotonList(fSpotPhotonList);

          ostringstream info1;
          info1 << "generating spot for: eye: " << fCurrentEyeId
                << ", telescope: " << fCurrentTelId
                << ", pixel: " << fCurrentPixId;
          INFO(info1);

          const double rDia = detTel.GetDiaphragmRadius();
          const double binWidth = detTel.GetCamera().GetFADCBinSize();
          const double normWavelength = detFD.GetReferenceLambda();

          const CoordinateSystemPtr& telCS = detTel.GetTelescopeCoordinateSystem();

          map<int, TraceD*> traces;
          for (int i = 0; i < fNSpotPhotons; ++i) {
            // Generating random point on the diaphragm
            const double diaPh = RandFlat::shoot(&fRandomEngine->GetEngine(), 0, 2*kPi);
            const double diaR = sqrt(RandFlat::shoot(&fRandomEngine->GetEngine(),
                                                     fRDiaMin*fRDiaMin, rDia*rDia));

            const double xDia = diaR*cos(diaPh);
            const double yDia = diaR*sin(diaPh);
            const Point pIn(xDia, yDia, 0, telCS);

            const double phi = detPixel.GetDirection().GetPhi(telCS);
            const double theta = detPixel.GetDirection().GetTheta(telCS);

            Vector nIn (1, theta, phi, telCS, Vector::kSpherical);
            nIn = -nIn;

            const double time = binWidth * RandFlat::shoot(&fRandomEngine->GetEngine(),
                                                           fMinBin, fMaxBin);
            const double weight = 1;
            utl::Photon photonIn(pIn, nIn, normWavelength, weight);
            photonIn.SetTime(TimeInterval(time));
            telSim.AddPhoton(photonIn);
          } // loop over raytraces
        } // end extra scope

        finishedOneSpot = true;
        return true;

      innerLoop:
        continue;

      } // end loop over pixels

      iPix = 1;
      FinishTelescope();

    } // end loop tels

  } // End loop over Eyes

  return false;
}


// -----------------------------------------------------------------------
void
SpotPhotonGenerator::FinishSpot()
{
  if (fNSpotPhotons <= 0) {
    ostringstream err;
    err << "Cannot continue with fNSpotPhotons=" << fNSpotPhotons;
    ERROR(err);
    return;
  }

  Detector& detector = Detector::GetInstance();
  const FDetector& detFD = detector.GetFDetector();
  const fdet::Pixel& detPixel = detFD.GetEye(fCurrentEyeId).GetTelescope(fCurrentTelId).GetPixel(fCurrentPixId);
  const CoordinateSystemPtr& pixelCS = detPixel.GetPixelCoordinateSystem();

  ostringstream spotname;
  spotname << "Eye_" << fCurrentEyeId << "_Tel_" << fCurrentTelId << "_Pix_" << fCurrentPixId;
  TH2D* const hSpot = new TH2D(spotname.str().c_str(), spotname.str().c_str(),
                               150, -2.75, 2.75, 150, -2.75, 2.75);
  hSpot->SetXTitle("X  [cm]");
  hSpot->SetYTitle("Y  [cm]");
  hSpot->SetZTitle("Photons");



  int nPhotons = 0;
  for (TelescopeSimulatorKG::TelescopeSimulator::PhotonListConstIterator iPh = fSpotPhotonList.begin();
       iPh != fSpotPhotonList.end(); ++iPh) {

    if (RandFlat::shoot(&fRandomEngine->GetEngine(), 0., 1.) > iPh->first.GetWeight())
      continue;

    ++nPhotons;
    hSpot->Fill(iPh->first.GetPosition().GetX(pixelCS)/cm,
                iPh->first.GetPosition().GetY(pixelCS)/cm);
  }

  const double efficiency = double(nPhotons) / fNSpotPhotons;
  const double efficiencyErr = sqrt((efficiency - efficiency*efficiency) / fNSpotPhotons);
  ostringstream info;
  info << " Efficiency is : " << setw(9) << setprecision(7) << efficiency
       << " +- " << efficiencyErr
       << " for " << fCurrentPixId;
  INFO(info);

  fEfficiencies[fCurrentPixId] = efficiency;
  fEfficienciesErr[fCurrentPixId] = efficiencyErr;

  gStyle->SetPalette(1);
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(1);
  TCanvas* const cSpot = new TCanvas("spot", "spot", 500, 500);
  cSpot->SetFillColor(0);
  hSpot->Draw("colz");
  spotname << ".eps";
  cSpot->Print(spotname.str().c_str());

  delete hSpot;
}


// -------------------------------------------------------------------
void
SpotPhotonGenerator::FinishTelescope()
{
  ostringstream filename;
  filename << "Efficiencies_Eye_" << fCurrentEyeId << "_Tel_" << fCurrentTelId << ".dat";
  ofstream file(filename.str().c_str());
  ostringstream info;
  info << " Pixel-Efficiencies for Eye: " << fCurrentEyeId << " Tel: " << fCurrentTelId << "\n";

  const int lineBreak = 10;

  double meanErr = 0;
  for (unsigned int i = 1; i <= fEfficiencies.size(); ++i) {
    meanErr += fEfficienciesErr[i];
    if ((i % lineBreak) == 1) {
      info << "       ";
      file << "       ";
    }
    info << setw(9) << setprecision(7) << fEfficiencies[i] << " ";
    file << setw(9) << setprecision(7) << fEfficiencies[i] << " ";
    if (((i + 1) % lineBreak) == 1) {
      info << endl;
      file << endl;
    }
  }

  if (fEfficiencies.size())
    meanErr /= fEfficiencies.size();

  info << "       <!-- Mean uncertainty on efficiencies is: " << meanErr << "  --> \n";
  file << "       <!-- Mean uncertainty on efficiencies is: " << meanErr << "  --> \n";
  info << "\n";
  file.close();
  INFO(info);

  fEfficiencies.clear();
  fEfficienciesErr.clear();
}


// -----------------------------------------------------------------------
VModule::ResultFlag
SpotPhotonGenerator::Finish()
{
  return eSuccess;
}


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