RdStationNoiseAnalyser.cc

Go to the documentation of this file.

#include "RdStationNoiseAnalyser.h"

#include <fwk/CentralConfig.h>
#include <fwk/RunController.h>
#include <fwk/VModule.h>

#include <evt/Event.h>
#include <revt/REvent.h>
#include <revt/Header.h>
#include <revt/Station.h>
#include <revt/Channel.h>
#include <utl/AugerUnits.h>
#include <utl/TimeStamp.h>
#include  <utl/ModifiedJulianDate.h>
#include <utl/UTCDateTime.h>
#include <utl/TimeInterval.h>
#include <utl/Branch.h>
#include<utl/SVector.h>
#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <rdet/Channel.h>
#include <rdet/Station.h>
#include <utl/AugerUnits.h>

#include <utl/ErrorLogger.h>
#include <string>
#include <complex>

#include <TH1F.h>
#include <TH2F.h>
#include <TProfile.h>
#include <TTimeStamp.h>
#include <TFile.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TGraph.h>

using namespace evt;
using namespace revt;
using namespace fwk;
using namespace utl;
using std::complex;
using std::vector;
using std::stringstream;
using std::ostringstream;


// Heres comes two function for converting time in LST
// No reason to keep them as class members
// Based on B. Revenu & G. Maurin Works

double
mod(const double d, const double periode)
{
  if (d >= 0)
    return d - floor(d/periode) * periode;
  else
    return d - ceil(d/periode) * periode + periode;
}


void
ct2lst(double longitude, double jd, double& lst)
{
  const double J2000 = 2451545.0;
  const double sidereal_coeff[4] = { 280.46061837, 360.98564736629, 0.000387933, 38710000 };

  double dt = jd - J2000;
  double dt0 = dt / 36525.;
  double theta = sidereal_coeff[0] + dt * sidereal_coeff[1] + dt0 * dt0 * (sidereal_coeff[2] - dt0 / sidereal_coeff[3]);        // is in degrees
  lst = mod(theta + longitude, 360.) * 0.066666666666666;
}


RdStationNoiseAnalyser::RdStationNoiseAnalyser() :
  fOutFile(NULL),
  fAllTrace(0),
  fLastEventTS(0, 0),
  fTimeMin(0),
  fTimeMax(0),
  fTSMin(0),
  fTSMax(0),
  fStationName("toto")
{
  //INFO("Creating RdStationNoiseAnalyser");
}


RdStationNoiseAnalyser::~RdStationNoiseAnalyser()
{
  delete fH_RateLST;
  delete fH_RMSLST1;
  delete fH_RMSLST2;
  delete fH_TimeSinceLastEvent;
  delete fH2_SpectrLST1;
  delete fH2_SpectrLST2;
  delete fH2_TimeSinceLastEventLST;
  delete fH_RateLocT;
  delete fH_RMSLocT2;
  delete fH_RMSLocT1;
  delete fH2_SpectrLocT1;
  delete fH2_SpectrLocT2;
  delete fH2_TimeSinceLastEventLocT;
  delete fH2_RMSLST2D1;
  delete fH2_RMSLST2D2;
  //delete fH2_RMSRate;
  delete fOutFile;
}


fwk::VModule::ResultFlag
RdStationNoiseAnalyser::Init()
{
  INFO("RdStationNoiseAnalyser::Init()");
  ostringstream info;
  Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdStationNoiseAnalyser");
  topBranch.GetChild("outputFile").GetData(fOutFileName);
  topBranch.GetChild("NBinTime").GetData(fNbinTime);
  topBranch.GetChild("NBinFreq").GetData(fNbinFreq);
  topBranch.GetChild("FrequencyMin").GetData(fFreqMin);
  topBranch.GetChild("FrequencyMax").GetData(fFreqMax);
  topBranch.GetChild("RMSMax").GetData(fRMSMax);
  topBranch.GetChild("PartOfTimeSeriesToUse").GetData(fPartOfTheTraceToUse);
  topBranch.GetChild("TimeMin").GetData(fTSMin);
  topBranch.GetChild("TimeMax").GetData(fTSMax);
  topBranch.GetChild("StationName").GetData(fStationName);
  float Ts_Min = float(fTSMin.GetGPSSecond() + kGPS_UTC_Offset);
  float Ts_Max = float(fTSMax.GetGPSSecond() + kGPS_UTC_Offset);
  ostringstream outinfo;
  outinfo << "\n-----------------------------------------------------------------------------\n" ;
  outinfo << int(fTSMin.GetGPSSecond()/*+kGPS_UTC_Offset*/) << " " << int(fTSMax.GetGPSSecond()/*+kGPS_UTC_Offset*/) << "\n";
  UTCDateTime totestmin(fTSMin.GetGPSSecond());
  outinfo << totestmin.GetYear() << "-"<<totestmin.GetMonth()<<"-"<<totestmin.GetDay() << "\n";
  outinfo << "-------------------------------------------------------------------------------\n";
  fNDay = TMath::Ceil((Ts_Max - Ts_Min)/(24*3600));
  INFO(outinfo.str().c_str());
  if (fNDay < 0) {
    ERROR("Error Number of day cannot be negative");
    return eFailure;
  }
  std::cout << fNDay << std::endl;
  //fOutFileName="/home/melissas/OffLine/MySeq/Noise.root";
  fH_RateLST = new TH1F("RateLST", "RateLST", fNbinTime, 0,24);
  fH_RateLocT = new TH1F("RateLocT", "RateLocT", fNbinTime, 0,24);
  fH_RateTime = new TH1F("RateTime", "RateTime", fNDay*fNbinTime, Ts_Min, Ts_Max);
  fH_RMSLST1 = new TH1F("RMSLST1", "RMSLST1", fNbinTime, 0,24);
  fH_RMSLST2 = new TH1F("RMSLST2", "RMSLST2", fNbinTime, 0,24);
  fH_RMSLocT1 = new TH1F("RMSLocT1", "RMSLocT1", fNbinTime, 0,24);
  fH_RMSLocT2 = new TH1F("RMSLocT2", "RMSLocT2", fNbinTime, 0,24);
  fH_TimeSinceLastEvent = new TH1F("TimeSinceLastEvent", "TimeSinceLastEvent", 200, 0, 4);
  fH2_TimeSinceLastEventLST = new TH2F("TimeSinceLastEventLST", "TimeSinceLastEventLST", fNbinTime, 0, 24, 800, 0, 4);
  fH2_TimeSinceLastEventLocT = new TH2F("TimeSinceLastEventLocT", "TimeSinceLastEventLocT", fNbinTime, 0, 24, 800, 0, 4);
  fH2_SpectrLST1 = new TH2F("SpectrLST1", "SpectrLST1", fNbinTime, 0,24, fNbinFreq, fFreqMin, fFreqMax);
  fH2_SpectrLST2 = new TH2F("SpectrLST2", "SpectrLST2", fNbinTime, 0,24, fNbinFreq, fFreqMin, fFreqMax);
  fH2_SpectrLocT1 = new TH2F("SpectrLocT1", "SpectrLocT1", fNbinTime, 0,24, fNbinFreq, fFreqMin, fFreqMax);
  fH2_SpectrLocT2 = new TH2F("SpectrLocT2", "SpectrLocT2", fNbinTime, 0,24, fNbinFreq, fFreqMin, fFreqMax);
  fH2_SpectrTime1 = new TH2F("SpectrTime1", "SpectrTime1", fNDay*fNbinTime, Ts_Min, Ts_Max, fNbinFreq, fFreqMin, fFreqMax);
  fH2_SpectrTime2 = new TH2F("SpectrTime2", "SpectrTime2", fNDay*fNbinTime, Ts_Min, Ts_Max, fNbinFreq, fFreqMin, fFreqMax);
  fH2_RMSLocT2D1 = new TH2F("RMSLocT2D1", "RMSLocT2D1", fNbinTime, 0, 24, fNDay, Ts_Min, Ts_Max);
  fH2_RMSLocT2D2 = new TH2F("RMSLocT2D2", "RMSLocT2D2", fNbinTime, 0, 24, fNDay, Ts_Min, Ts_Max);
  fH2_RateLocT2D = new TH2F("RateLocT2D", "RateLocT2D", fNbinTime, 0, 24, fNDay, Ts_Min, Ts_Max);
  fH2_RMSLST2D1 = new TH2F("RMSLST2D1", "RMSLST2D1", fNbinTime, 0, 24, 100, 0, fRMSMax);
  fH2_RMSLST2D2 = new TH2F("RMSLST2D2", "RMSLST2D2", fNbinTime, 0, 24, 100, 0, fRMSMax);
  if (fPartOfTheTraceToUse > 1.0)
    return eFailure;

  info << " Data Will be saved to \n" << fOutFileName
       << " Number Of Time Bin (0h-24h)" << fNbinTime << "\n"
          " Number Of Freq Bin (" << fFreqMin << " MHz - " << fFreqMax << " MHz) : " << fNbinFreq << "\n"
          " Will Use the last " << fPartOfTheTraceToUse*100 << "% of the time series";
  INFO(info);
  return eSuccess;
}


fwk::VModule::ResultFlag
RdStationNoiseAnalyser::Run(evt::Event& theevent)
{
  const rdet::RDetector& rDetector =
    det::Detector::GetInstance().GetRDetector();
  revt::REvent& therev = theevent.GetREvent();
  vector<double> v_RMS1;
  vector<double> v_RMS2;
  v_RMS1.clear();
  v_RMS2.clear();
  double lst = 0;
  double loctime = 0;
  if (fTimeMin == 0 && fTimeMax == 0) {
    fTimeMin = therev.GetHeader().GetTime().GetGPSSecond();
    fTimeMax = therev.GetHeader().GetTime().GetGPSSecond();
  } else {
    if (fTimeMin > therev.GetHeader().GetTime().GetGPSSecond())
      fTimeMin = therev.GetHeader().GetTime().GetGPSSecond();
    if (fTimeMax < therev.GetHeader().GetTime().GetGPSSecond())
      fTimeMax = therev.GetHeader().GetTime().GetGPSSecond();
  }
  TimeInterval interval = therev.GetHeader().GetTime() - fLastEventTS;
  UTCDateTime utcdate(therev.GetHeader().GetTime());

  int EventSecond = therev.GetHeader().GetTime().GetGPSSecond() + kGPS_UTC_Offset;
  fLastEventTS = therev.GetHeader().GetTime();
  ct2lst(-69, utl:: ModifiedJulianDate(utcdate), lst);
  loctime = double(utcdate.GetHour()) + double(utcdate.GetMinute())/60.;
  fH_RateLST->Fill(lst);
  fH_RateLocT->Fill(loctime);
  fH_RateTime->Fill(EventSecond);
  if (interval.GetNanoSecond() > 0 && interval.GetSecond() < 4) {
    double delay = double(interval.GetSecond()) + double(interval.GetNanoSecond())/1e9;
    fH_TimeSinceLastEvent->Fill(delay);
    fH2_TimeSinceLastEventLST->Fill(lst, delay);
    fH2_TimeSinceLastEventLocT->Fill(loctime, delay);
  }
  try {
    Station& stat=therev.GetStationByName(fStationName);
    if (!stat.HasChannel(0) && !stat.HasChannel(1)) {
      WARNING("Channel 0 or 1 not present \n");
    }
    ConvertChannelToStation(stat, rDetector);
    //Creating a virtual station with the physical amplitude
//    StationFFTDataContainer virtualstationdata;
//    StationFrequencySpectrum& virtualstationspectrum =
//      virtualstationdata.GetStationFrequencySpectrum();
//    Vector3C V3CZero;
//    Vector3C efield;
//    V3CZero = complex<double>(0.0), complex<double>(0.0), complex<double>(0.0);
//    for (int ichan=0;ichan<=1;ichan++) {

  //    }
    const StationTimeSeries& stseries = stat.GetStationTimeSeries();
    StationTimeSeries theseries;
    StationFFTDataContainer stationData;
    //float normfactorsp=1;
    int binmin = floor((1 - fPartOfTheTraceToUse)*stseries.GetSize());
    for (unsigned int i = binmin; i < stseries.GetSize(); ++i) {
      //double Value[3] = { stseries[i][0], stseries[i][1], stseries[i][2] };
      SVector<3, double> Value;
      Value = stseries[i][0], stseries[i][1], stseries[i][2];
      v_RMS1.push_back(Value[0]/(micro*volt/meter));
      v_RMS2.push_back(Value[1]/(micro*volt/meter));
      theseries.PushBack(Value);
    }
    theseries.SetBinning(stseries.GetBinning());
    stationData.GetTimeSeries() = theseries;
    stationData.SetNyquistZone(stat.GetNyquistZone());
    const revt::StationFrequencySpectrum& stspectr = stationData.GetFrequencySpectrum();
    float normfactorsp = 0;
//    float normfactorsp = abs(float(stspectr.GetSize()) / float(stat.GetFrequencyOfBin(0) - stat.GetFrequencyOfBin(stspectr.GetSize()-1)))*micro*volt/meter/megahertz;

   //  normfactorsp = abs(float(stat.GetFrequencyOfBin(0) - stat.GetFrequencyOfBin(stspectr.GetSize()-1)))/float(stspectr.GetSize());
    normfactorsp = micro*volt/meter/megahertz;
//    normfactorsp/=(micro*volt/meter);
    for (unsigned int i = 0; i < stspectr.GetSize(); ++i) {
      fH2_SpectrLST1->Fill(lst, stationData.GetFrequencyOfBin(i)/megahertz, abs(stspectr[i][0])/normfactorsp) ;
      fH2_SpectrLST2->Fill(lst, stationData.GetFrequencyOfBin(i)/megahertz, abs(stspectr[i][1])/normfactorsp);
      fH2_SpectrLocT1->Fill(loctime, stationData.GetFrequencyOfBin(i)/megahertz, abs(stspectr[i][0])/normfactorsp);
      fH2_SpectrLocT2->Fill(loctime, stationData.GetFrequencyOfBin(i)/megahertz, abs(stspectr[i][1])/normfactorsp);
      fH2_SpectrTime1->Fill(EventSecond, stationData.GetFrequencyOfBin(i)/megahertz, abs(stspectr[i][0])/normfactorsp);
      fH2_SpectrTime2->Fill(EventSecond, stationData.GetFrequencyOfBin(i)/megahertz, abs(stspectr[i][1])/normfactorsp);
    }
    //cout << therev.GetHeader().GetTime().GetGPSSecond() << " // " << TimeDay.GetGPSSecond() <<endl;
    double Rms1 = TMath::RMS(v_RMS1.size(), &v_RMS1.front());
    double Rms2 = TMath::RMS(v_RMS2.size(), &v_RMS1.front());
    fH2_RMSLocT2D1->Fill(loctime, EventSecond-3600*loctime, Rms1);
    fH2_RMSLocT2D2->Fill(loctime, EventSecond-3600*loctime, Rms2);
    fH2_RateLocT2D->Fill(loctime, EventSecond-3600*loctime);
    fH_RMSLST1->Fill(lst,Rms1);
    fH_RMSLST2->Fill(lst,Rms2);
    fH2_RMSLST2D1->Fill(lst, Rms1);
    fH2_RMSLST2D2->Fill(lst, Rms2);
    fH_RMSLocT1->Fill(loctime, Rms1);
    fH_RMSLocT2->Fill(loctime, Rms2);
  } catch (utl::NonExistentComponentException& e) {
    ostringstream warn;
    warn.str("");
    warn << "Station " << fStationName << "does not exist for event " << therev.GetHeader().GetId();
    WARNING(warn.str());
  }
  return eSuccess;
}


fwk::VModule::ResultFlag
RdStationNoiseAnalyser::Finish()
{
  INFO("RdStationNoiseAnalyser::Finish");
  INFO("Writing outfile");
  float TotTime = float(fTimeMax - fTimeMin);
  fOutFile = new TFile(fOutFileName.c_str(),"RECREATE");
//  fH_RateLST->Scale(1./TotTime);
  fH_RateLST->Write(); // Can still be used for normalisation (even /s ) but what with empty bin ?
  fH_RMSLST1->Write();
  fH_RMSLST2->Write();
  fH_TimeSinceLastEvent->Write();
  fH2_SpectrLST1->Write();
  fH2_SpectrLST2->Write();
  fH2_TimeSinceLastEventLST->Write();
  fH_RateLocT->Scale(1./(TotTime));
  fH_RateLocT->Write();
  fH_RMSLocT1->Write();
  fH_RMSLocT2->Write();
  fH2_SpectrLocT1->Write();
  fH2_SpectrLocT2->Write();
  fH2_SpectrTime1->Write();
  fH2_SpectrTime2->Write();
  fH_RateTime->Write();
  fH2_RateLocT2D->Write();
  fH2_RMSLocT2D1->Write();
  fH2_RMSLocT2D2->Write();
  fH2_RMSLST2D1->Write();
  fH2_RMSLST2D2->Write();
  fH2_RMSLocT2D2->Write();
  fH2_TimeSinceLastEventLST->Write();
  TProfile* Prof_TimeSinceLastEventLST = fH2_TimeSinceLastEventLST->ProfileX("TSinceLastEventProfLST");
  Prof_TimeSinceLastEventLST->Write();
  delete Prof_TimeSinceLastEventLST;
  fH2_TimeSinceLastEventLocT->Write();
  TProfile* Prof_TimeSinceLastEventLocT = fH2_TimeSinceLastEventLocT->ProfileX("TSinceLastEventProfLocT");
  Prof_TimeSinceLastEventLocT->Write();
  delete Prof_TimeSinceLastEventLocT;
  fOutFile->Close();
  delete fOutFile;
  fOutFile = 0;
  ostringstream info;
  info << "Total data Time " << TotTime << " (s) ";
  INFO(info.str());
  return eSuccess;
}
void RdStationNoiseAnalyser::ConvertChannelToStation(revt::Station& stat,const rdet::RDetector& rDet) {
      // StationTimeSeries& newstationtimeseries=stat.GetStationTimeSeries();  // unused variable
      StationFrequencySpectrum& stspectr=stat.GetStationFrequencySpectrum();
      bool firststationchannel = true;
      Vector3C V3CZero;
      V3CZero= complex<double>(0.0),complex<double>(0.0),complex<double>( 0.0);
      for (Station::ChannelIterator cIt = stat.ChannelsBegin();
           cIt != stat.ChannelsEnd(); ++cIt) {
        Channel& channel = *cIt;
  if (channel.GetId()>2) continue ;
/*      fMessage.str("");
        fMessage << "Copying channel "
                 << channel.GetId()
                 << " of station "
                 << channel.GetStationId()
                 << " to the StationTimeSeries.";
        INFO(fMessage.str()); */

        // Work on the trace of this channel
      //  const revt::ChannelTimeSeries& channeltrace =
        //  channel.GetConstChannelTimeSeries();  // unused
  const revt::ChannelFrequencySpectrum& channelspectr=
    channel.GetConstChannelFrequencySpectrum();
//  TGraph gr_chanresponse( channelspectr.GetSize());
        // when working on the first channel, fill station timeseries with zeroes and set some station-specific quantities
        if (firststationchannel) {
          firststationchannel = false;
          // overwrite the station time series with an empty one of the correct size
           stspectr  = StationFrequencySpectrum(channelspectr.GetSize(),channelspectr.GetBinning(),V3CZero);
          // set the Nyquist zone of the station, does no check if Nyquist zones of Channels are the same!
          stat.SetNyquistZone(channel.GetNyquistZone());
        }

        // check if newstationtimeseries.GetSize() is correct
        if (stspectr.GetSize() != channelspectr.GetSize()) {
          stringstream fMessage;
          fMessage << "Number of samples in spectrum of "
                   << channel.GetId()
                   << " belonging to station "
                   << channel.GetStationId()
                   << " is different than the number of samples in other channels of that station.";
          WARNING(fMessage.str());
        }

        // now copy over the values, assume Channel 1 is x, Channel 2 is y, Channel 3 is z
        for (ChannelFrequencySpectrum::SizeType i=0; i<channelspectr.GetSize(); ++i) {
    double channelfrequency=channel.GetFrequencyOfBin(i);
//          newstationtimeseries[i][channel.GetId()-1] = channeltrace[i];
    double response=  rDet.GetStation(stat.GetId()).GetChannel(channel.GetId()).GetIntegratedEffectiveAntennaHeight(channelfrequency);
//    gr_chanresponse.SetPoint(i,channelfrequency/megahertz,response);
    if (response==0) {
      ostringstream warn;
      warn.str("");
      warn << " Detector response is zero !! \n";
      warn << " Frequency is " << channelfrequency/megahertz ;
      WARNING(warn.str());
      continue;
      }
    stspectr[i][channel.GetId()-1]=channelspectr[i]/response;
//    cout << "VALUES >>>>>> " << channelspectr[i] << " " << response << channelspectr[i]/response << endl;
        }
    //  ostringstream title;
   //   title.str("");
   //   title << "Response_Station"<<rDet.GetStation(stat.GetId()).GetName()  <<"_Chan_"<<channel.GetId();
   //   gr_chanresponse.SetTitle(title.str().c_str());
   //   TCanvas canrep(title.str().c_str(),title.str().c_str(),640,480);
    //  gr_chanresponse.Draw("AL");
     // canrep.Print((title.str()+".ps").c_str());
      }

  }