RdChannelGalacticConstantsGenerator.cc

Go to the documentation of this file.

#include "RdChannelGalacticConstantsGenerator.h"
#include <utl/ErrorLogger.h>

#include <fwk/CentralConfig.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <revt/StationRecData.h>
#include <revt/StationGPSData.h>

#include <rdet/RDetector.h>
#include <rdet/Channel.h>

#include <utl/Trace.h>
#include <utl/TraceAlgorithm.h>
#include <utl/Reader.h>
#include <utl/config.h>
#include <utl/AugerUnits.h>
#include <utl/SaveCurrentTDirectory.h>

#include <TFile.h>
#include <TTree.h>
#include <TH2D.h>
#include <TProfile2D.h>

#include <boost/scoped_ptr.hpp>

#include <limits>


using namespace std;
using namespace utl;
using namespace fwk;
using namespace revt;
using namespace rdet;


namespace RdChannelGalacticConstantsGenerator {


  VModule::ResultFlag
  RdChannelGalacticConstantsGenerator::Init()
  {
    INFO("RdChannelGalacticConstantsGenerator::Init()");

    Branch topBranch =
      CentralConfig::GetInstance()->GetTopBranch("RdChannelGalacticConstantsGenerator");

    topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
    topBranch.GetChild("IgnoreStations").GetData(fIgnoreStations);
    topBranch.GetChild("MaxIntensity").GetData(fMaxIntensity);
    topBranch.GetChild("OutputFileName").GetData(fOutputFileName);
    topBranch.GetChild("SaveROOTFile").GetData(fSaveROOTFile);
    topBranch.GetChild("FileNameROOTFile").GetData(fFileNameROOTFile);

    fOutfile.open(fOutputFileName.c_str());

    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelGalacticConstantsGenerator::Run(evt::Event& event)
  {
    INFO("RdChannelGalacticConstantsGenerator::Run()");

    REvent& rEvent = event.GetREvent();

    for (auto& station : rEvent.CandidateStationsRange()) {
      const int stationId = station.GetId() - 100;
      const unsigned long gpsSec = station.GetGPSData().GetSecond();

      if (gpsSec < fLowestGPSSec)
        fLowestGPSSec = gpsSec;
      if (gpsSec > fHighestGPSSec)
        fHighestGPSSec = gpsSec;

      for (auto& channel : station.ChannelsRange()) {

        if (!channel.IsActive())
          continue;

        const int channelId = channel.GetId();
        const rdet::RDetector& rDetector = det::Detector::GetInstance().GetRDetector();
        const rdet::Channel& detChannel = rDetector.GetStation(station.GetId()).GetChannel(channel.GetId());
        const double lowerDesignFreq = detChannel.GetDesignLowerFreq();
        const double upperDesignFreq = detChannel.GetDesignUpperFreq();

        const revt::ChannelFrequencySpectrum& spectrum = channel.GetChannelFrequencySpectrum();

        const int lowerBin = int(lowerDesignFreq/spectrum.GetBinning());
        const int upperBin = int(upperDesignFreq/spectrum.GetBinning());
        const int nSpecBins = upperBin - lowerBin;

        const double lowerFreq = channel.GetFrequencyOfBin(lowerBin + 1) / utl::MHz;
        const double upperFreq = channel.GetFrequencyOfBin(upperBin + 1) / utl::MHz;

        const unsigned int fSpectrumSize = spectrum.GetSize();

        // calculate lst hour of current event
        const double lst = GPSSecToLSTHourAuger(gpsSec);

        // some hard-coded binning
        const int stationNumber = 200; // not very likely to increase in the forseeable future
        const int lstBins = 144; // 10-minute bins, should be enough

        // first event, create TProfiles to fill data
        if (fFirstEvent) {
          fLSTSpecNS = new TProfile3D("LSTSpecNS", "LSTSpecNS",
                                      stationNumber, 0, stationNumber,
                                      lstBins, 0, 24,
                                      nSpecBins, lowerFreq, upperFreq);
          fLSTSpecEW = new TProfile3D("LSTSpecEW", "LSTSpecEW",
                                      stationNumber, 0, stationNumber,
                                      lstBins, 0, 24,
                                      nSpecBins, lowerFreq, upperFreq);
          fFirstEvent = false;
        }

        double intensity = 0;
        int intensityEntries = 0;

        if (fLSTSpecNS->GetNbinsZ() != nSpecBins) {
          ostringstream warn;
          warn << "Number of spectral bins (" << nSpecBins
               << ") does not correspond to profile bins (" << fLSTSpecNS->GetNbinsZ()
               << "). Input data probably has multiple trace lengths.";
          WARNING(warn);
        }

        // check if predefined normalized intensity cut is not crossed
        for (ChannelFrequencySpectrum::SizeType i = 0; i < fSpectrumSize; ++i) {
          const double currentFreq = channel.GetFrequencyOfBin(i) / utl::MHz;
          if (lowerFreq <= currentFreq && currentFreq < upperFreq) {
            intensity += abs(spectrum[i])*abs(spectrum[i]);
            ++intensityEntries;
          }
        }
        if (sqrt(intensity) / intensityEntries > fMaxIntensity) {
          ostringstream warn;
          warn << "Average intensity above cut, skipping this trace from station "
               << stationId + 100 << " channel " << channelId;
          WARNING(warn);
          continue;
        }

        // fill values in TProfiles
        for (ChannelFrequencySpectrum::SizeType i = 0; i < fSpectrumSize; ++i) {
          if (channelId == 1) {
            const double currentFreq = channel.GetFrequencyOfBin(i) / utl::MHz;
            if (lowerFreq <= currentFreq && currentFreq < upperFreq) {
              fLSTSpecNS->Fill(stationId, lst, currentFreq + 0.5*spectrum.GetBinning() / utl::MHz,
                               abs(spectrum[i]));
            }
          }
          if (channelId == 2) {
            const double currentFreq = channel.GetFrequencyOfBin(i) / utl::MHz;
            if (lowerFreq <= currentFreq && currentFreq < upperFreq) {
              fLSTSpecEW->Fill(stationId, lst, currentFreq + 0.5*spectrum.GetBinning() / utl::MHz,
                               abs(spectrum[i]));
            }
          }
        }
      } // channel loop
    } // station loop
    return eSuccess;
  }


  VModule::ResultFlag
  RdChannelGalacticConstantsGenerator::Finish()
  {
    INFO("RdChannelGalacticConstantsGenerator::Finish()");

    // If option is selected, write profiles to root tree and file for later analysis
    if (fSaveROOTFile) {
      SaveCurrentTDirectory save;
      const boost::scoped_ptr<TFile> file(new TFile(fFileNameROOTFile.c_str(), "RECREATE"));
      TTree* const tree = new TTree("fTree", "data");
      tree->Branch("pLSTSpecNS", "TProfile3D", &fLSTSpecNS);
      tree->Branch("pLSTSpecEW", "TProfile3D", &fLSTSpecEW);
      tree->Branch("fLowestGPSSec", &fLowestGPSSec);
      tree->Branch("fHighestGPSSec", &fHighestGPSSec);
      tree->Fill();
      tree->Write();
      file->Close();
    }

    // This part of the module will process the profiles and calculate the gain factors //

    // project the profile to count the number of recorded stations, and see which stations correspond to stations to be ignored
    vector<int> vActiveStations;
    int nIgnoreCorrection = 0;
    const boost::scoped_ptr<TH1D> hLSTSpecNSprojX(fLSTSpecNS->ProjectionX("hLSTSpecNSprojX"));
    std::ostringstream info;
    info << "Active station IDs for which calibration constants will be calculated: ";
    for (int jj = 1; jj <= hLSTSpecNSprojX->GetNbinsX(); ++jj) {
      if (hLSTSpecNSprojX->GetBinContent(jj)) {
        vActiveStations.push_back(jj - 1);
        info << jj+99 << ' ';
        for (unsigned int iignore = 0; iignore < fIgnoreStations.size(); ++iignore) {
          if (jj-1 == fIgnoreStations[iignore]-100) {
            ++nIgnoreCorrection;
            info << "(ignored for average) ";
          }
        }
      }
    }
    INFO(info);
    info.str("");
    const int nActiveStations = vActiveStations.size();

    // project the profile to count the number of filled LST bins
    const int totalLSTbins = fLSTSpecNS->GetNbinsY();
    const boost::scoped_ptr<TH2D> hLSTSpecNSprojLST((TH2D*)fLSTSpecNS->Project3D("xy"));
    hLSTSpecNSprojLST->SetName("hLSTSpecNSprojLST");
    for (int jj = 0; jj < nActiveStations; ++jj) {
      int lstBins = 0;
      for (int ybin = 1; ybin <= hLSTSpecNSprojLST->GetNbinsY(); ++ybin) {
        const int bin = hLSTSpecNSprojLST->GetBin(ybin, vActiveStations[jj]+1);
        if (hLSTSpecNSprojLST->GetBinContent(bin))
          ++lstBins;
      }
      const double lstFraction= double(lstBins) / totalLSTbins;
      if (lstFraction < 1) {
        std::ostringstream warn;
        warn << "Data from station " << vActiveStations[jj]+100
             << " does not fully cover LST range (fraction = " << lstFraction
             << "). This will result in a wrong calibration!";
        WARNING(warn);
      }
    }

    // set band stopped bins to zero and make bins empty in profiles
    const boost::scoped_ptr<TProfile3D> pLSTSpecNSFiltered((TProfile3D*)fLSTSpecNS->Clone());
    delete fLSTSpecNS;
    fLSTSpecNS = nullptr;
    const boost::scoped_ptr<TProfile3D> pLSTSpecEWFiltered((TProfile3D*)fLSTSpecEW->Clone());
    delete fLSTSpecEW;
    fLSTSpecEW = nullptr;
    pLSTSpecNSFiltered->SetName("LSTSpecNSFiltered");
    pLSTSpecEWFiltered->SetName("LSTSpecEWFiltered");
    for (int iz = 1; iz <= pLSTSpecNSFiltered->GetNbinsZ(); ++iz) {
      for (int ix = 1; ix <= pLSTSpecNSFiltered->GetNbinsX(); ++ix) {
        for (int iy = 1; iy <= pLSTSpecNSFiltered->GetNbinsY(); ++iy) {
          const int currentbin = pLSTSpecNSFiltered->GetBin(ix, iy, iz);
          if (pLSTSpecNSFiltered->GetBinContent(currentbin) < 1e-9) {
            pLSTSpecNSFiltered->SetBinContent(currentbin, 0);
            pLSTSpecNSFiltered->SetBinEntries(currentbin, 0);
          }
        }
      }
    }
    for (int iz = 1; iz <= pLSTSpecEWFiltered->GetNbinsZ(); ++iz) {
      for (int ix = 1; ix <= pLSTSpecEWFiltered->GetNbinsX(); ++ix) {
        for (int iy = 1; iy <= pLSTSpecEWFiltered->GetNbinsY(); ++iy) {
          const int currentbin = pLSTSpecEWFiltered->GetBin(ix, iy, iz);
          if (pLSTSpecEWFiltered->GetBinContent(currentbin) < 1e-9) {
            pLSTSpecEWFiltered->SetBinContent(currentbin, 0);
            pLSTSpecEWFiltered->SetBinEntries(currentbin, 0);
          }
        }
      }
    }

    // clone and remove ignored stations to calculate average
    const boost::scoped_ptr<TProfile3D> pLSTSpecNSGoodStations((TProfile3D*)pLSTSpecNSFiltered->Clone());
    const boost::scoped_ptr<TProfile3D> pLSTSpecEWGoodStations((TProfile3D*)pLSTSpecEWFiltered->Clone());
    pLSTSpecNSGoodStations->SetName("LSTSpecNSGoodStations");
    pLSTSpecEWGoodStations->SetName("LSTSpecEWGoodStations");
    for (unsigned int iignore = 0; iignore < fIgnoreStations.size(); ++iignore) {
      for (int iz = 1; iz <= pLSTSpecNSGoodStations->GetNbinsZ(); ++iz) {
        for (int iy = 1; iy <= pLSTSpecNSGoodStations->GetNbinsY(); ++iy) {
          int zerobin = pLSTSpecNSGoodStations->GetBin(fIgnoreStations[iignore]-99, iy, iz);
          pLSTSpecNSGoodStations->SetBinContent(zerobin, 0);
          pLSTSpecNSGoodStations->SetBinEntries(zerobin, 0);
          zerobin = pLSTSpecEWGoodStations->GetBin(fIgnoreStations[iignore]-99, iy, iz);
          pLSTSpecEWGoodStations->SetBinContent(zerobin, 0);
          pLSTSpecEWGoodStations->SetBinEntries(zerobin, 0);
        }
      }
    }

    // create total average LST spectra by projecting along the x-axis (station axis)
    const boost::scoped_ptr<TH2D> lstSpecNSFilteredzy((TH2D*)pLSTSpecNSGoodStations->Project3D("zy"));
    lstSpecNSFilteredzy->SetName("LSTSpecNSFiltered_zy");
    lstSpecNSFilteredzy->Scale(1. / (nActiveStations - nIgnoreCorrection));
    const boost::scoped_ptr<TH2D> lstSpecEWFilteredzy((TH2D*)pLSTSpecEWGoodStations->Project3D("zy"));
    lstSpecEWFilteredzy->SetName("LSTSpecEWFiltered_zy");
    lstSpecEWFilteredzy->Scale(1. / (nActiveStations - nIgnoreCorrection));

    info << "The following averaged calibration constants were calculated (NS, EW):";

    ostringstream str1;
    for (int jj = 0; jj < nActiveStations; ++jj) {

      // create complete and filtered LST spectra of individual stations
      pLSTSpecNSFiltered->GetXaxis()->SetRange(vActiveStations[jj]+1, vActiveStations[jj]+1); // in tprofile, add one to get correct bin
      pLSTSpecEWFiltered->GetXaxis()->SetRange(vActiveStations[jj]+1, vActiveStations[jj]+1); // set range to limit profile to one particular station
      str1.str("");
      str1 << "zy_stat" << vActiveStations[jj] << " ns", vActiveStations[jj];
      // project the profile to get the spectrum of the one station in its range
      const boost::scoped_ptr<TProfile2D> pLSTSpecNSFilteredDivided((TProfile2D*)pLSTSpecNSFiltered->Project3D(str1.str().c_str()));
      str1.str("");
      str1 << "zy_stat" << vActiveStations[jj] << " ew";
      const boost::scoped_ptr<TProfile2D> pLSTSpecEWFilteredDivided((TProfile2D*)pLSTSpecEWFiltered->Project3D(str1.str().c_str()));
      pLSTSpecNSFilteredDivided->Divide(lstSpecNSFilteredzy.get());// divide the projection by the total average lst spectrum profile (this/average)
      pLSTSpecEWFilteredDivided->Divide(lstSpecEWFilteredzy.get());

      // now calculate divided spectra for single stations, these histograms are the frequency-dependent strengths as measured
      const boost::scoped_ptr<TH1D> hCalibrationConstantsNS(pLSTSpecNSFilteredDivided->ProjectionY()); // project onto LST-axis
      const boost::scoped_ptr<TH1D> hCalibrationConstantsEW(pLSTSpecEWFilteredDivided->ProjectionY());
      hCalibrationConstantsNS->Scale(1. / totalLSTbins); // normalise to adjust for number of LST bins
      hCalibrationConstantsEW->Scale(1. / totalLSTbins);

      // output frequency-dependent calibration constants of the stations to file
      for (int ibin = 1; ibin <= hCalibrationConstantsNS->GetNbinsX(); ++ibin) {
        fOutfile << vActiveStations[jj] << "   1   " << hCalibrationConstantsNS->GetBinCenter(ibin)
                 << "   " << hCalibrationConstantsNS->GetBinContent(ibin) << endl;
      }
      for (int ibin = 1; ibin <= hCalibrationConstantsEW->GetNbinsX(); ++ibin) {
        fOutfile << vActiveStations[jj] << "   2   " << hCalibrationConstantsEW->GetBinCenter(ibin)
                 << "   " << hCalibrationConstantsEW->GetBinContent(ibin) << endl;
      }

      // calculate and output average calibration constants for quick overview
      info << "\n  station ID " << vActiveStations[jj]+100 << ": ";
      double averageConstantsNS = 0;
      int constentries = 0;
      for (int ibin = 0; ibin < hCalibrationConstantsNS->GetNbinsX(); ++ibin) {
        if (hCalibrationConstantsNS->GetBinContent(ibin)) {
          averageConstantsNS += hCalibrationConstantsNS->GetBinContent(ibin);
          ++constentries;
        }
      }
      averageConstantsNS /= constentries;
      double averageConstantsEW = 0;
      constentries = 0;
      for (int ibin = 0; ibin < hCalibrationConstantsEW->GetNbinsX(); ++ibin) {
        if (hCalibrationConstantsEW->GetBinContent(ibin)) {
          averageConstantsEW += hCalibrationConstantsEW->GetBinContent(ibin);
          ++constentries;
        }
      }
      averageConstantsEW /= constentries;
      info << '(' << averageConstantsNS << ", " << averageConstantsEW << ")";

    }

    INFO(info);

    fOutfile.close();

    return eSuccess;
  }


#define J2000 2451545.0
#define Dtohr 0.066666666666666
#define Hrtod 15

  // function to calculate LST
  double
  RdChannelGalacticConstantsGenerator::GPSSecToLSTHourAuger(const int gps)
  {
     // FIXME DV: this code is horrible. who added this crap? this has to be removed ASAP.
    const time_t utc = gps + 315964819 + 16; // GPS UTC offset, Leap Seconds have to be updated
    const tm* const t = gmtime(&utc);
    int year = t->tm_year + 1900;
    int month = t->tm_mon + 1;
    const int day = t->tm_mday;
    const int hour = t->tm_hour;
    const int min = t->tm_min;
    const int sec = t->tm_sec;

//#warning can you use the utl::ModifiedJulianDate here?
    //Date to Julian Date
    if (month <= 2) {
      --year;
      month += 12;
    }
    const int a = int(floor(year / 100.));
    const int b = 2 - a + a/4;
    double jd = floor(365.25*(year + 4716)) + floor(30.6001*(month + 1)) + day + b - 1524.5;
    const double fd = double(60*60*hour + 60*min + sec) / 86400;
    jd = jd + fd;

    static double sidereal_coeff[4] = { 280.46061837, 360.98564736629, 0.000387933, 38710000 };
    const double dt = jd - J2000;
    const double dt0 = dt / 36525.;
    const double theta =
      sidereal_coeff[0] + dt * sidereal_coeff[1] +
      dt0 * dt0 * (sidereal_coeff[2] - dt0 / sidereal_coeff[3]); // is in degrees
    // original J. Meeus: longitude is positive westwards  *lst = mod(theta-longitude,360.) * kSTC::Dtohr;
    // with the new convention for longitudes (positive eastwards) :
    const double lst = fmod(theta + (360 - 69.6), 360.) * Dtohr;
    return lst;
  }

}