TelescopeSimulatorLX/TelescopeSimulator.cc

Go to the documentation of this file.

#include <utl/config.h>

#include "TelescopeSimulator.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 <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 <CLHEP/Random/Randomize.h>

#include <boost/tuple/tuple.hpp>

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

#include <G4RunManager.hh>
#include <G4UImanager.hh>
#include <G4VUserPhysicsList.hh>
#include <G4Timer.hh>
#include <G4ProcessTable.hh>
//#include <G4VisExecutive.hh>

#include "FDsimG4DetectorConstruction.hh"
#include "FDsimG4PrimaryGeneratorAction.hh"
#include "FDsimG4PhysicsList.hh"
#include "FDsimG4EventAction.hh"
#include "FDsimG4StoreOpticalHit.hh"
#include "FDsimG4VisExecutive.hh"
#include "FDsimG4XMLManager.hh"
#include "FDsimG4Defs.hh"
#include "FDsimG4Write2ROOT.hh"

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

using CLHEP::RandFlat;


VModule::ResultFlag
TelescopeSimulator::Init()
{
  fXMLManager = FDsimG4XMLManager();
  fXMLManager.Init("TelescopeSimulatorLX");

  fDrumMode = false;

  INFO("Geant4 telescope ray-tracing (beta-version)");

  Branch topB = CentralConfig::GetInstance()->GetTopBranch("TelescopeSimulatorLX");
  if (!topB) {
    ERROR("No XML configuration found for TelescopeSimulatorLX!");
    return eFailure;
  }

  // ---------------------  **OPTIONAL**  ---------------------
  fVerbosityLevel = 0;

  if (topB.GetChild("verbosityLevel"))
    topB.GetChild("verbosityLevel").GetData(fVerbosityLevel);

  fG4MacroFile = "";
  fG4Standalone = false;
  fWeightFactor = 1.0;
  fROOTFile  = "FDsimG4.root";

  Branch G4optionsB = topB.GetChild("options");
  if (G4optionsB) {
    string instr;

    if (G4optionsB.GetChild("ROOTFile"))
      G4optionsB.GetChild("ROOTFile").GetData(fROOTFile);

    if (G4optionsB.GetChild("G4MacroFile"))
      G4optionsB.GetChild("G4MacroFile").GetData(fG4MacroFile);

    if (G4optionsB.GetChild("G4Standalone")) {
      G4optionsB.GetChild("G4Standalone").GetData(instr);
      fG4Standalone = (instr == "yes");
    }

    if (G4optionsB.GetChild("weightFactor"))
      G4optionsB.GetChild("weightFactor").GetData(fWeightFactor);
  }

  if (!fG4Standalone)
    fWeightFactor = 1.0;
  if (fG4Standalone) {
    INFO("TelescopeSimulatorLX will be run in standalone mode.");
    if (fG4MacroFile == "") {
      ERROR ("Standalone mode was selected but no macro file given.");
      return eFailure;
    }
  }

  if (fG4MacroFile != "")
    INFO("Reading Geant4 macro file.");

  fWrite2ROOT = nullptr;

#ifdef __TELESCOPE_SIMULATOR_LX_USE_ROOT__
  if (!fROOTFile.empty())
    fWrite2ROOT = new FDsimG4Write2ROOT(fROOTFile);
#endif

  INFO("Initializing FDsimG4");

  fG4RunManager = new G4RunManager();

  FDsimG4DetectorConstruction* const myDetector =
    new FDsimG4DetectorConstruction("CO1", fXMLManager);
  myDetector->SetVerbosityLevel(fVerbosityLevel);
  fG4RunManager->SetUserInitialization(myDetector);
  fG4RunManager->SetUserInitialization(new FDsimG4PhysicsList);

  FDsimG4PrimaryGeneratorAction* fG4PrimGenAct = new FDsimG4PrimaryGeneratorAction();
  fG4RunManager->SetUserAction(fG4PrimGenAct);

  fG4PrimGenAct->SetUseGeneralParticleSource(fG4Standalone);

  if (fWrite2ROOT)
    fG4EventAction = new FDsimG4EventAction(fWrite2ROOT);
  else
    fG4EventAction = new FDsimG4EventAction();

  fG4RunManager->SetUserAction(fG4EventAction);

  FDsimG4VisExecutive* const visManager = new FDsimG4VisExecutive();
  visManager->Initialize();

  // Initialize G4 kernel
  fG4RunManager->Initialize();

  // Get the Pointer to the UI Manager
  fG4UImanager = G4UImanager::GetUIpointer();

  /* --------------------------------------------------- the model
     config signature
  */
//#warning This is purposely hardcoded to 'TelescopeSimulatorKG' in order to allow the approximate calibration of the raytracing. TelescopeSimulatorLX is not to be used in any real simulation production nor physical analysis like this.
  fGeneralConfigSignature = "TelescopeSimulatorKG";

  return eSuccess;
}


VModule::ResultFlag
TelescopeSimulator::Run(evt::Event& event)
{
  if (fG4Standalone && !event.HasFEvent())
    event.MakeFEvent();

  if (fG4Standalone) {

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

    const FDetector& detFD = detector.GetFDetector();

    const TimeStamp& simTime = detector.GetTime();

    // set the FD event times
    fEvent.GetHeader().SetTime(simTime);
    event.GetHeader().SetTime(simTime);
    event.GetHeader().SetId("");

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

      const int eyeId = itEye->GetId();

      if (!event.GetFEvent().HasEye(eyeId))
        event.GetFEvent().MakeEye(eyeId);

      fevt::Eye& eyeEvent = event.GetFEvent().GetEye(eyeId);

      for (Eye::TelescopeIterator itTel = itEye->TelescopesBegin();
           itTel != itEye->TelescopesEnd(); ++itTel) {

        const int telId = itTel->GetId();

        if (!eyeEvent.HasTelescope(telId))
          eyeEvent.MakeTelescope(telId);
        if (!eyeEvent.GetTelescope(telId).HasSimData())
          eyeEvent.GetTelescope(telId).MakeSimData();

        fevt::Telescope& telData = eyeEvent.GetTelescope(telId);

        fevt::TelescopeSimData& telSim = telData.GetSimData();
        telSim.SetPhotonsStartTime(simTime);

        const unsigned int telTraceNBins = 1000;
        telSim.SetNumberOfPhotonBins(telTraceNBins);

      }

    }

  }

  //

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

  if (!event.HasSimShower()) {
    /* In drum 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;
    if (fWeightFactor != 1) {
      ERROR("Photon weight-factor was not 1. Reseting!");
      fWeightFactor = 1;
    }
  }

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

  //

  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;

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

      Reset();

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

      const CoordinateSystemPtr& telCS = detTel.GetTelescopeCoordinateSystem();

#ifdef __TELESCOPE_SIMULATOR_LX_USE_ROOT__
      const fdet::Eye& detEye = detFD.GetEye(*iEye);
      const CoordinateSystemPtr& eyeCS = detEye.GetEyeCoordinateSystem();
      if (fWrite2ROOT) {
        fWrite2ROOT->SetTelCS(telCS);
        fWrite2ROOT->SetEyeCS(eyeCS);
      }
#endif

      if (fVerbosityLevel > 0) {
        ostringstream info;
        info << "Raytracing eye=" << eyeId << ", "
                "telescope=" << telId << ", "
                "rDia=" << setw(6) << detTel.GetDiaphragmRadius();
        INFO(info);
      }

      if (!fG4MacroFile.empty()) {
        const std::string command = "/control/execute ";
        fG4UImanager->ApplyCommand(command + fG4MacroFile);
      }

      if (!fG4Standalone) {

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

        const int nPhotons = telSim.GetNPhotons();
        if (!nPhotons)
          continue;

        if (fVerbosityLevel > 0) {
          ostringstream info;
          info << "Raytracing nPhotons=" << nPhotons;
          INFO(info);
        }

        FDsimG4StoreOpticalHit* g4Photon = nullptr;

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

        for (fevt::TelescopeSimData::PhotonIterator iPhoton = telSim.PhotonsBegin();
             iPhoton != telSim.PhotonsEnd(); ++iPhoton) {

          const double wavelength = iPhoton->GetWavelength();

          if (wavelength < minWl || wavelength > maxWl)
            continue;

          utl::Photon photonIn = *iPhoton;

          g4Photon = new FDsimG4StoreOpticalHit();

          const auto wl = photonIn.GetWavelength();
          g4Photon->SetEnergy(utl::kPlanck * utl::kSpeedOfLight / wl / utl::eV);
          g4Photon->SetWavelength(wl / utl::nanometer);
          g4Photon->SetTime(photonIn.GetTime().GetInterval() / utl::ns);

          const auto& pos = photonIn.GetPosition();
          g4Photon->SetPosX(pos.GetX(telCS) / utl::m);
          g4Photon->SetPosY(pos.GetY(telCS) / utl::m);
          g4Photon->SetPosZ(pos.GetZ(telCS) / utl::m);

          const auto& dir = photonIn.GetDirection();
          g4Photon->SetDirX(dir.GetX(telCS));
          g4Photon->SetDirY(dir.GetY(telCS));
          g4Photon->SetDirZ(dir.GetZ(telCS));

          g4Photon->SetWeight(photonIn.GetWeight());

          fG4Photons.push_back(g4Photon);

        }

        const int nPhot = fG4Photons.size();

        ostringstream info;
        info << "Processing " << nPhot << " photons...";
        INFO(info);

        fG4PhotonsIter = fG4Photons.begin();

        fG4RunManager->BeamOn(nPhot);

      }

      FillTraces(iEye, iTel);

    }

  }

  return eSuccess;
}


VModule::ResultFlag
TelescopeSimulator::Finish()
{
  if (fG4RunManager)
    delete fG4RunManager ;

  G4VVisManager * visManager  = G4VVisManager::GetConcreteInstance();
  if (visManager)
    delete visManager;

  cout << "Vis manager deleted." << endl;

  if (fG4EventAction)
    fG4EventAction = NULL;

#ifdef __TELESCOPE_SIMULATOR_LX_USE_ROOT__
  if (fWrite2ROOT){
    fWrite2ROOT->SaveAndClose();
    delete fWrite2ROOT;
    fWrite2ROOT = NULL;
  }
#endif

  return eSuccess;
}


FDsimG4StoreOpticalHit*
TelescopeSimulator::GetG4Photon()
{
  FDsimG4StoreOpticalHit* const g4Phot = *fG4PhotonsIter;
  ++fG4PhotonsIter;
  return g4Phot;
}


void
TelescopeSimulator::SetG4PixelHit(const FDsimG4StoreOpticalHit& g4PixelHit)
{
  fG4PixelHits.push_back(g4PixelHit);
}


void
TelescopeSimulator::Reset()
{
  // CBT Reset here vectors
  for (const auto p : fG4Photons)
    delete p;

  fG4Photons.clear();
  fG4PixelHits.clear();
}


void
TelescopeSimulator::FillTraces(const fevt::FEvent::EyeIterator& iEye,
                               const fevt::Eye::TelescopeIterator& iTel)
{
  Detector& detector = Detector::GetInstance();
  const FDetector& detFD = detector.GetFDetector();

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

  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();

  const double normWavelength = detFD.GetReferenceLambda();

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

  map<int, TraceD*> traces;

  FDsimG4StoreOpticalHit pixelHit;

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

  for (const auto& pixelHit : fG4PixelHits) {

    const unsigned int pxlId = pixelHit.GetPMTid();
    const double wavelength  = pixelHit.GetWavelength() * utl::nanometer;
    const double time = pixelHit.GetTime() * utl::ns;
    const double weight = pixelHit.GetWeight();

    if (fVerbosityLevel > 1)
      pixelHit.Print();

    if (wavelength < minWl || wavelength > maxWl)
      continue;

    const int col = (pxlId - 1) / 22 + 1;
    const int row = pxlId - 22*(col - 1) ;

    if (pxlId < detTel.GetFirstPixelId() || pxlId > detTel.GetLastPixelId()) {
      ostringstream err;
      err << " Pixel Id=" << pxlId << " does not exists on camera!"
          << " col=" << col << " row=" << row;
      ERROR(err);
      continue;
    }

    const fdet::Pixel& detPix = detTel.GetPixel(pxlId);
    const TabulatedFunction& qEff = detPix.GetQEfficiency();

    if (wavelength < qEff[0].X() || wavelength > qEff[qEff.GetNPoints()-1].X())
      continue;

    // !! this factor is to achive consistency to reconstruction !!
    double efficiencyCorrection = 0;
    if (!fDrumMode) {
      const double modelRelEff = detTel.GetModelRelativeEfficiency(wavelength);
      if (modelRelEff)
        efficiencyCorrection = telMeasEfficiency.Y(wavelength) / modelRelEff;
    }

    const int bin = int(time / telTraceBinWidth);

    if (fVerbosityLevel > 1) {
      cerr << " bin=" << setw(4) << bin
           << " time=" << setw(6) << time/utl::ns
           << " weight=" << setw(4) << weight
           << " wl=" << setw(4) << wavelength/utl::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;
    }

    if (bin < 0 || bin >= int(telTraceNBins)) {

      if (fVerbosityLevel > 2) {
        ostringstream err;
        err << "Bin out of range. bin=" << bin
            << " weight=" << weight << "telTraceNBins= " << telTraceNBins;
        ERROR(err);
      }

      continue;
    }

    // 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] +=
      (fDrumMode ? 1 : (weight * fWeightFactor
                        * qEff.Y(wavelength) / qEff.Y(normWavelength)
                        * efficiencyCorrection));

  }

  if (fVerbosityLevel > 2) {
    // Output pixels with signal
    for (auto iTrace = traces.begin(); iTrace != traces.end(); ++iTrace) {

      const int pxlId = iTrace->first;

      const int col = (pxlId - 1) / 22 + 1; // (col-1)*22 + row;
      const int row = (pxlId - 1) % 22 + 1;

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

      double signal = 0;
      for (unsigned int i = 0; i < telTraceNBins; ++i) {

        if (fVerbosityLevel > 3 && (*iTrace->second)[i]) {
          cerr << " i: " << i << " val: " << (*iTrace->second)[i]
               << " |";
        }
        signal += (*iTrace->second)[i];
      }

      cerr << " total signal " << signal << endl;

    }
  }
}