DrumPhotonGenerator.cc

Go to the documentation of this file.

#include "DrumPhotonGenerator.h"

#include <utl/config.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/SaveCurrentTDirectory.h>

#include <fwk/CentralConfig.h>
#include <fwk/CoordinateSystemRegistry.h>
#include <fwk/RandomEngineRegistry.h>
#include <fwk/GITGlobalRevision.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 <evt/LaserData.h> // Modified by D'Urso & Valore 02/03/06

#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 <TDirectory.h>
#include <TFile.h>
#include <TTree.h>

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

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

using CLHEP::RandFlat;


DrumPhotonGenerator::DrumPhotonGenerator()
{ }


DrumPhotonGenerator::~DrumPhotonGenerator()
{ }


VModule::ResultFlag
DrumPhotonGenerator::Init()
{
  CentralConfig* cc = CentralConfig::GetInstance();
  Branch topB = cc->GetTopBranch("DrumPhotonGenerator");

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

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

  // Initialize the drum simulation mode

  topB.GetChild("minAccuracy").GetData(fMinAccuracy);
  topB.GetChild("NDrumPhotons").GetData(fNDrumPhotons);
  topB.GetChild("nCos").GetData(fDrumNCos);
  topB.GetChild("MinThetaPhotons").GetData(fThetaMinDrumPhotons);
  topB.GetChild("MaxThetaPhotons").GetData(fThetaMaxDrumPhotons);
  topB.GetChild("MinPhiPhotons").GetData(fPhiMinDrumPhotons);
  topB.GetChild("MaxPhiPhotons").GetData(fPhiMaxDrumPhotons);
  topB.GetChild("MinRDiaPhotons").GetData(fRDiaMin);
  topB.GetChild("MaxBin").GetData(fMaxBin);
  topB.GetChild("MinBin").GetData(fMinBin);
  topB.GetChild("DrumDataOut").GetData(fDrumDataOutName);

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

  // info output
  ostringstream info;
  info << " Version: "
       << GetVersionInfo(VModule::eRevisionNumber) << "\n"
          " Parameters:\n"
          "         min accuracy: " << fMinAccuracy << "\n"
          "    number of photons: " << fNDrumPhotons << "\n"
          "                 nCos: " << fDrumNCos;
  if (fDrumNCos == 0)
    info << " (isotropic)\n";
  if (fDrumNCos == 1)
    info << " (lambertian)\n";
  info << "            min theta: " << setw(3) << fThetaMinDrumPhotons/deg << "deg\n"
          "            max theta: " << setw(3) << fThetaMaxDrumPhotons/deg << "deg\n"
          "              min phi: " << setw(3) << fPhiMinDrumPhotons/deg << "deg\n"
          "              max phi: " << setw(3) << fPhiMaxDrumPhotons/deg << "deg\n"
          "           min radius: " << setw(3) << fRDiaMin/m << "m\n"
          "              min bin: " << fMinBin << "\n"
          "              max bin: " << fMaxBin << "\n"
          "     drum-data output: " << fDrumDataOutName;
  INFO(info);

  fInit = false;
  fStatus = eGeneratePhotons;

  return eSuccess;
}


VModule::ResultFlag
DrumPhotonGenerator::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();

  if (!fInit) {
    INFO("Creating FEvent");
    fInit = true;

    const TimeStamp& simTime = detector.GetTime();

    // set the FD event times
    fEvent.GetHeader().SetTime(simTime);
    event.GetHeader().SetTime(simTime);
    ostringstream idStr;
    //idStr << "DrumSimulation_" << 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");
        }

        const int index = Index(eyeId, telId);
        fTotalGeneratedPhotonsPerTelescope[index] = 0;

        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 (CalculateCalibrationConstants(event))
    return eSuccess;
  else if (DoDrum(event))
    return eSuccess;
  return eBreakLoop;
}


// -----------------------------------------------------------------------
bool
DrumPhotonGenerator::DoDrum(evt::Event& event)
{
  fevt::FEvent& fEvent = event.GetFEvent();

  bool accuracyNotYetReached = false;

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

    for (fevt::Eye::TelescopeIterator iTel = iEye->TelescopesBegin(fevt::ComponentSelector::eInDAQ);
         iTel != iEye->TelescopesEnd(fevt::ComponentSelector::eInDAQ); ++iTel) {
      const int telId = iTel->GetId();
      const int index = Index(eyeId, telId);

      fevt::TelescopeSimData& telSim = iTel->GetSimData();
      telSim.ClearPhotons();

      ostringstream info1;
      info1 << "calibrating eye: " << eyeId << ", telescope: " << telId
            << ", nTotalPhotons: " << fTotalGeneratedPhotonsPerTelescope[index] ;
      INFO(info1);

      map<int, CalibResult> constants = CalibrateTelescope(*iTel);

      // check for current precision (-> photon statistics at the pixels)
      double meanAccuracy = 100.;
      double minAccuracy = 100.;
      double maxAccuracy = 100.;
      double meanCal = 100.;
      double minCal = 100.;
      double maxCal = 100.;
      bool firstAccuracy = true;

      int countPixels = 0;
      for (map<int, CalibResult>::const_iterator iCal = constants.begin();
           iCal != constants.end(); ++iCal) {

        if (iCal->second.calib <= 0)
          continue;

        ++countPixels;
        const double accuracy = iCal->second.calibErr / iCal->second.calib;

        if (firstAccuracy) {
          firstAccuracy = false;
          meanAccuracy = accuracy;
          minAccuracy = accuracy;
          maxAccuracy = accuracy;
          meanCal = iCal->second.calib;
          minCal = iCal->second.calib;
          maxCal = iCal->second.calib;
        } else {
          meanAccuracy += accuracy;
          if (accuracy < minAccuracy)
            minAccuracy = accuracy;
          if (accuracy > maxAccuracy)
            maxAccuracy = accuracy;
          meanCal += iCal->second.calib;
          if (iCal->second.calib < minCal)
            minCal = iCal->second.calib;
          if (iCal->second.calib > maxCal)
            maxCal = iCal->second.calib;
        }
      }

      if (countPixels) {
        meanAccuracy /= countPixels;
        meanCal /= countPixels;

        ostringstream info2, info3;
        info2 << "-> accuracy (min/max/mean): " << minAccuracy << " / " << maxAccuracy << " / " << meanAccuracy;
        info3 << "-> calib (min/max/mean): " << minCal << " / " << maxCal << " / " << meanCal;
        INFO(info2);
        INFO(info3);
      }

      if (maxAccuracy > fMinAccuracy) {
        telSim.ClearConfigSignature();
        GenerateDrumPhotons(*iTel);
        accuracyNotYetReached = true;
      } else {
        ostringstream info2;
        info2 << "===> Accuracy in this telescope is reached. Switching it off now. ";
        INFO(info2);

        // this telescope is done -> take it out of DAQ
        iTel->SetStatus(fevt::ComponentSelector::eExists);
      }

    } // end loop tels

  } // End loop over Eyes

  if (accuracyNotYetReached)
    return true;

  INFO("Drum photon precision achieved.");
  fStatus = eCalibrate;

  return false; // calibration need no more photons
}


void
DrumPhotonGenerator::GenerateDrumPhotons(fevt::Telescope& tel)
{
  ostringstream info;
  info << "Generating " << fNDrumPhotons << " photons from drum";
  INFO(info);

  const unsigned int telId = tel.GetId();
  const int index = Index(tel.GetEyeId(), telId);

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

  const fdet::Telescope& detTel = detFD.GetTelescope(tel);
  const double rDia = detTel.GetDiaphragmRadius();
  const double binWidth = detTel.GetCamera().GetFADCBinSize();
  const double normWavelength = detFD.GetReferenceLambda();

  fevt::TelescopeSimData& telSim = tel.GetSimData();

  const CoordinateSystemPtr& telCS = detTel.GetTelescopeCoordinateSystem();

  fTotalGeneratedPhotonsPerTelescope[index] += fNDrumPhotons;

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

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

    // Generating random direction on the diaphragm
    double phi = RandFlat::shoot(&fRandomEngine->GetEngine(),
                                 fPhiMinDrumPhotons, fPhiMaxDrumPhotons);
    double ncos = fDrumNCos; // th goes like cos^NCOS * d(cos)
    double theta_1 = std::pow(cos(fThetaMinDrumPhotons), ncos+1);
    double theta_2 = std::pow(cos(fThetaMaxDrumPhotons), ncos+1);
    // th goes like cos^NCOS * d(cos)
    double tmp = (theta_2-theta_1) * RandFlat::shoot(&fRandomEngine->GetEngine(),
                                                     0.0, 1.0) + theta_1;

    double theta;
    if (ncos == 0)      theta = acos(tmp);
    else if (ncos == 1) theta = acos(sqrt(tmp));
    else                theta = acos(exp(log(tmp) / (1.+ncos) ) );

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

    double time = binWidth * RandFlat::shoot(&fRandomEngine->GetEngine(), fMinBin, fMaxBin);
    const double weight = 1.0;
    utl::Photon photonIn(pIn, nIn, normWavelength, weight);
    photonIn.SetTime(TimeInterval(time));
    telSim.AddPhoton(photonIn);

  } // loop over raytraces
}


bool
DrumPhotonGenerator::CalculateCalibrationConstants(evt::Event& event)
{
  if (fStatus == eGeneratePhotons)
    return false;
  if (fStatus == eCalibrated) {
    ERROR(" CANNOT RE-CALIBRATE CALIBRATED TELESCOPE !!!!!!!!!!!!!!!!!!!!! ");
    return false;
  }
  fStatus = eCalibrated;

  INFO("Calibrate telescopes ....");

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

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

  ofstream calibXml("FSimulationCalibConfig.generated.xml.in");
  calibXml << "<?xml version=\"1.0\" encoding=\"iso-8859-1\"?>\n\n"
              "<FSimulationCalibConfig\n"
              "  xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"\n"
              "  xsi:noNamespaceSchemaLocation='@XMLSCHEMALOCATION@/FSimulationCalibConfig.xsd'>\n\n"
              "  <!-- these are used for the flatfielding -->\n\n";

  const SaveCurrentTDirectory save;
  TFile* const drumDataFile = TFile::Open(fDrumDataOutName.c_str(), "UPDATE");
  drumDataFile->cd();
  TTree* const runInfo = new TTree("runInfo", "runInfo");

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

    for (fevt::Eye::TelescopeIterator iTel = iEye->TelescopesBegin(fevt::ComponentSelector::eExists);
         iTel != iEye->TelescopesEnd(fevt::ComponentSelector::eExists); ++iTel) {
      const int telId = iTel->GetId();
      const int index = Index(eyeId, telId);

      const fevt::TelescopeSimData& telSim = iTel->GetSimData();
      const string configList = telSim.GetConfigSignatureStr();
      const string md5 = telSim.GetConfigSignature();

      calibXml << "  <!-- \n"
                  "    This set of calibration constants was generated using Offline revision "
               << fwk::GITGlobalRevision::GetInstance().GetId() << "\n"
                  "    with the following set of configuration parameters:\n"
                  "\"" << configList << "\"\n"
                  "    The following Md5-hex-digest encodes this list:\n"
                  "    -->\n"
                  "  <EndToEnd signature=\"" << md5 << "\"> \n"
                  "    <calibration>" << endl;

      drumDataFile->cd();
      ostringstream treeName;
      treeName << "drumData_" << index;
      TTree* const drumData = new TTree(treeName.str().c_str(), treeName.str().c_str());
      int drumData_id;
      double drumData_nGen, drumData_nPix, drumData_calib, drumData_calibErr, drumData_elecFact;
      drumData->Branch("id", &drumData_id, "id/I");
      drumData->Branch("nGen", &drumData_nGen, "nGen/D");
      drumData->Branch("nPix", &drumData_nPix, "nPix/D");
      drumData->Branch("calib", &drumData_calib, "calib/D");
      drumData->Branch("calibErr", &drumData_calibErr, "calibErr/D");
      drumData->Branch("electronicsFactor", &drumData_elecFact, "electronicsFactor/D");

      map<int, CalibResult> constants = CalibrateTelescope(*iTel);

      const fdet::Telescope& detTel = detFD.GetTelescope(*iTel);
      const double diaArea = detTel.GetDiaphragmArea();

      // const double solidAngle = 2. * M_PI * (1. - cos(fThetaMaxDrumPhotons));
      const double meanPixelSolid = 5.94e-4 * sr;

      // flux of drum vertical to axis (for lambertian)
      const double flux_drum_perp =
        fTotalGeneratedPhotonsPerTelescope[index] / (diaArea * M_PI * (1 - pow(cos(fThetaMaxDrumPhotons), 2)));

      ostringstream info1, info2;
      info1 << " Eye: " << eyeId << ", Telescope: " << telId
            << ", Has lens: " << (detTel.HasCorrectorRing() ? "yes" : "no")
            << ", Signature: " << md5 << "\n"
               "total number of ph.: " << fTotalGeneratedPhotonsPerTelescope[index] << " photons\n"
               "    perp. drum flux: " << flux_drum_perp * cm2 * sr <<  "1/(cm^2 deg^2)" << "\n"
               " mean pixel photons: " << flux_drum_perp * diaArea * meanPixelSolid << " photons ";
      INFO(info1);
      info2 << '\n';

      // copy first to local string not to take address of a temporary
      const string tree_rev = fwk::GITGlobalRevision::GetInstance().GetId();
      if (!runInfo->GetBranch("index")) {
        runInfo->Branch("index", const_cast<int*>(&index), "index/I");
        runInfo->Branch("nPhotonsTot", &fTotalGeneratedPhotonsPerTelescope[index], "nPhotonsTot/i");
        runInfo->Branch("md5", const_cast<char*>(md5.c_str()), "md5/C");
        runInfo->Branch("offline_rev", const_cast<char*>(tree_rev.c_str()), "offline_rev/C");
        runInfo->Branch("config", const_cast<char*>(configList.c_str()), "config/C");
      } else {
        runInfo->SetBranchAddress("index", const_cast<int*>(&index));
        runInfo->SetBranchAddress("nPhotonsTot", &fTotalGeneratedPhotonsPerTelescope[index]);
        runInfo->SetBranchAddress("md5", const_cast<char*>(md5.c_str()));
        runInfo->SetBranchAddress("offline_rev", const_cast<char*>(tree_rev.c_str()));
        runInfo->SetBranchAddress("config", const_cast<char*>(configList.c_str()));
      }

      runInfo->Fill();

      for (unsigned int iPix = detTel.GetFirstPixelId(); iPix <= detTel.GetLastPixelId(); ++iPix) {

        const int pixelId = iPix;

        const double calib = constants[pixelId].calib;
        const double calibErr = constants[pixelId].calibErr;

        info2 << " pixel id=" << setw(4) << pixelId
              << " s=" << setw(10) << fixed << calib
              << " +- " << setw(10) << calibErr  << " [ph/ADC]"
                 "\n";

        calibXml << "      " << setprecision(16) << setw(18) << calib
                 << "  <!-- +- " << setw(18) << calibErr << " -->\n";

        drumData_id = pixelId;
        drumData_nGen = constants[pixelId].nPhotGen;
        drumData_nPix = constants[pixelId].nPhotPixel;
        drumData_calib = constants[pixelId].calib;
        drumData_calibErr = constants[pixelId].calibErr;
        drumData_elecFact = constants[pixelId].electronicsFactor;
        drumData->Fill();

      } // end loop pixels

      INFO(info2);

      calibXml << "    </calibration>\n"
                  "  </EndToEnd>\n\n";

      // drumData->Write();

    } // end loop telescopes
  } // end loop eyes

  calibXml << "</FSimulationCalibConfig>\n";

  //runInfo->Write();

  drumDataFile->Write();
  drumDataFile->Close();

  INFO("Telescope(s) calibrated! Output in file \"FSimulationCalibConfig.generated.xml.in\".");

  return true;
}


map<int, DrumPhotonGenerator::CalibResult>
DrumPhotonGenerator::CalibrateTelescope(fevt::Telescope& tel)
{
  const int telId = tel.GetId();
  const int eyeId = tel.GetEyeId();
  const int index = Index(eyeId, telId);

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

  const fdet::Telescope& detTel = detFD.GetTelescope(tel);
  const double diaArea = detTel.GetDiaphragmArea();

  // flux of drum vertical to axis (for lambertian)
  const double flux_drum_perp = fTotalGeneratedPhotonsPerTelescope[index] /
    (diaArea * M_PI * (1 - pow(cos(fThetaMaxDrumPhotons), 2)));

  map<int, CalibResult> constants;
  if (fTotalGeneratedPhotonsPerTelescope[index] == 0)
    return constants;

  for (unsigned int iPix = detTel.GetFirstPixelId();
       iPix <= detTel.GetLastPixelId(); ++iPix) {

    const fevt::Pixel& pixel = tel.GetPixel(iPix);
    const fdet::Pixel& detPixel = detFD.GetPixel(pixel);
    const fdet::Channel& detChannel = detFD.GetChannel(pixel); // mapped channel

    const double pixelSolid = detPixel.GetSolidAngle();
    const double qEff = detPixel.GetQEfficiency().Y(normWavelength);
    const double gain = detChannel.GetElectronicsGain();

    const fevt::PixelSimData& pixelSimData = pixel.GetSimData();
    const TraceD& photonTrace = pixelSimData.GetPhotonTrace(fevt::FdConstants::eTotal);

    double pixelPhotons = 0;
    for (unsigned int i = 0; i < photonTrace.GetSize(); ++i)
      pixelPhotons += photonTrace[i];

    // this is really: nGenPhotons/cos(theta)
    const double nGenPhotons = flux_drum_perp * diaArea * pixelSolid;

    const double epsilon = pixelPhotons / nGenPhotons;
    const double epsilonErr = sqrt((epsilon - epsilon*epsilon) / nGenPhotons);

    // const double pixelSignal = pixelPhotons * qEff * gain; // add electronics [ph -> ADC]
    // const double calib = flux_drum_perp * diaArea * pixelSolid / pixelSignal;

    const double calib = 1 / (epsilon * qEff * gain);
    const double calibErr = epsilonErr / (epsilon * epsilon * qEff * gain);

    constants[iPix].calib = calib;
    constants[iPix].calibErr = calibErr;
    constants[iPix].nPhotGen = nGenPhotons;
    constants[iPix].nPhotPixel = pixelPhotons;
    constants[iPix].electronicsFactor = qEff * gain;
  }
  return constants;
}


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


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