RdChannelNoiseASCIIImporter.cc

Go to the documentation of this file.

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

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

#include <det/Detector.h>
#include <rdet/RDetector.h>
#include <rdet/Station.h>
#include <rdet/Channel.h>

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

#include <vector>
#include <map>
#include <algorithm>
#include <fwk/CentralConfig.h>

#include <time.h>
#include <sstream>

#include <boost/tokenizer.hpp>
#include <boost/algorithm/string/replace.hpp>
using namespace std;

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

typedef boost::tokenizer<boost::char_separator<char> > mytok;

namespace RdChannelNoiseASCIIImporter {

RdChannelNoiseASCIIImporter::RdChannelNoiseASCIIImporter() :
    fNumberOfFiles(0),
    fInfoLevel(0),
    fPathToAsciitrace(""),
    fUseAlwaysTheSameTrace(false),
    fUseRandomSeed(false),
    fStopAtLastNoiseTrace(true)
{
}

RdChannelNoiseASCIIImporter::~RdChannelNoiseASCIIImporter()
{
}

VModule::ResultFlag RdChannelNoiseASCIIImporter::Init()
{

  // Initialize your module here. This method
  // is called once at the beginning of the run.
  // The eSuccess flag indicates the method ended
  // successfully.  For other possible return types,
  // see the VModule documentation.

  INFO("RdChannelNoiseASCIIImporter::Init()");

  Branch topBranch = CentralConfig::GetInstance()->GetTopBranch("RdChannelNoiseASCIIImporter");
  topBranch.GetChild("NoiseFilePath").GetData(fPathToAsciitrace);
  topBranch.GetChild("InfoLevel").GetData(fInfoLevel);
  topBranch.GetChild("NumberOfTraces").GetData(fNumberOfFiles);
  topBranch.GetChild("UseAlwaysTheSameTrace").GetData(fUseAlwaysTheSameTrace);
  topBranch.GetChild("UseRandomSeed").GetData(fUseRandomSeed);
  topBranch.GetChild("StopAtLastNoiseTrace").GetData(fStopAtLastNoiseTrace);
  //Add a / at the end of the path if missing
  if (fPathToAsciitrace.compare(fPathToAsciitrace.size() - 1, fPathToAsciitrace.size(), "/") != 0)
    fPathToAsciitrace.append("/");
  if (fUseRandomSeed) {
    std::srand(time(NULL));
  }
  //Randomize the available traces
  for (int i = 1; i <= fNumberOfFiles; i++) {
    v_id.push_back(i);
  }
  std::random_shuffle(v_id.begin(), v_id.end());
  viter = v_id.begin();
  return eSuccess;
}

VModule::ResultFlag RdChannelNoiseASCIIImporter::Run(evt::Event& event)
{
  INFO("RdChannelNoiseASCIIImporter::Run()");
  if (viter == v_id.end()) {
    if (!fStopAtLastNoiseTrace) {
      WARNING("end of noise list reached, traces will be reshuffled and used again");
      std::random_shuffle(v_id.begin(), v_id.end());
      viter = v_id.begin();
    } else {
      cerr << "End of the Noise List" << endl;
      return eContinueLoop;
    }
  }
  vector<int> v_stationid;
  ostringstream fname;
  fname.str("");
  fname << fPathToAsciitrace << "Traces_" << AddZero(*viter) << *viter << ".dat";
  INFO(fname.str());
  ifstream inf(fname.str().c_str());
  while (!inf.good()) {
    INFO("file does not exist, moving to next file");
    ++viter;
    if (viter == v_id.end()) {
      if (!fStopAtLastNoiseTrace) {
        WARNING("end of noise list reached, traces will be reshuffled and used again");
        std::random_shuffle(v_id.begin(), v_id.end());
        viter = v_id.begin();
      } else {
        cerr << "End of the Noise List" << endl;
        return eContinueLoop;
      }
    }
    fname.str("");
    fname << fPathToAsciitrace << "Traces_" << AddZero(*viter) << *viter << ".dat";
    INFO(fname.str());
    inf.open(fname.str().c_str());
  }
  if (!fUseAlwaysTheSameTrace) {
    ++viter;  // We go here to the next value of viter
  }
  string zeline;
  boost::char_separator<char> sep(" ");
  map<string, vector<int> > m_traces;
  m_traces.clear();
  int currentStationId  = 0;
  while (getline(inf, zeline)) {
    vector<int> v_values;
    vector<string> v_tokens;
    // Just in case of clear the vector
    v_values.clear();
    v_tokens.clear();
    mytok tok(zeline, sep);  //Split the line into tok;
    for (mytok::const_iterator titer = tok.begin(); titer != tok.end(); ++titer) {  //Put all tokens in a vector;
      v_tokens.push_back(*titer);
    }
    for (unsigned int i = 3; i < v_tokens.size(); i++) {
      int val = atoi(v_tokens[i].c_str());
      v_values.push_back(val);
    }
    string name = GetStringFromStatChan(atoi(v_tokens[0].c_str()), atoi(v_tokens[1].c_str()));
    if(atoi(v_tokens[0].c_str()) != currentStationId) {
      v_stationid.push_back(atoi(v_tokens[0].c_str()));
      currentStationId = atoi(v_tokens[0].c_str());
    }
    RemovePedestal(v_values);
    m_traces[name] = v_values;
  }
  inf.close();
  random_shuffle(v_stationid.begin(), v_stationid.end());
  vector<int>::iterator randomstationiterator = v_stationid.begin();

  // Check if there are events at all
  if (!event.HasREvent()) {
    WARNING("RdChannelNoiseASCIIImporter::No radio event found!");
    return eContinueLoop;
  } else {
    det::Detector& det = det::Detector::GetInstance();
    ostringstream out;
    REvent& rEvent = event.GetREvent();
    // loop through stations and for every station through every channel
    bool randomNoiseStationUsed = false;
    // reshuffle random station iterator before each event to prevent selecting the same noise trace twice in one event
    for (REvent::StationIterator sIt = rEvent.StationsBegin(); sIt != rEvent.StationsEnd(); ++sIt) {
      Station& station = *sIt;
      if(randomNoiseStationUsed) {
        std::advance(randomstationiterator, 1); // increase random station iterator by 1 to get a different station next time
        randomNoiseStationUsed = false;
      }
      for (Station::ChannelIterator cIt = station.ChannelsBegin(); cIt != station.ChannelsEnd();
          ++cIt) {
        Channel& channel = *cIt;
        revt::ChannelTimeSeries& timeSeries = channel.GetChannelTimeSeries();

        string name = GetStringFromStatChan(station.GetId(), channel.GetId());
        map<string, vector<int> >::const_iterator miter = m_traces.find(name);
        if (miter == m_traces.end()) {
          if (randomstationiterator == v_stationid.end()) {
            std::random_shuffle(v_stationid.begin(), v_stationid.end());
            randomstationiterator = v_stationid.begin();
//#warning TH: this randomization is dangerous because the sample rate can be different for different stations and this is not checked for!
          }
          randomNoiseStationUsed = true;
          out.str("");
          out << " Unable to Find Station_Channel " << name;
          name = GetStringFromStatChan(*randomstationiterator, channel.GetId());
          out << " -> Will use a station " << name << " instead";
          if (fInfoLevel > eFinal) {
            INFO(out.str());
          }
          miter = m_traces.find(name);
          if (miter == m_traces.end()) {
            if (fInfoLevel > eFinal) {
              out.str("");
              out << " Damned we are unlucky station " << name << " is not found I give up. Setting channel as inactive";
              WARNING(out.str());
            }
            channel.SetNotActive();
            continue;
          }
        }

        // get number of adc counts from RDetector
        short bitDepth = det.GetRDetector().GetStation(station.GetId()).GetChannel(channel.GetId())
            .GetBitDepth();
        short numberOfADCCounts = pow(2., bitDepth) - 1;

        if (miter->second.size() != timeSeries.GetSize()) {
          out.str("");
          out << " Station_Channel " << name << "\n";
          out << "Channel size " << timeSeries.GetSize() << " \n";
          out << "Noise size " << miter->second.size() << "\n";
          out << " Channel " << name
              << " Is set as inactive, please check if this is coherent with what you are doing";
          WARNING(out.str().c_str());
          channel.SetNotActive();
          continue;
        }
        for (unsigned int i = 0; i < timeSeries.GetSize(); ++i) {
          int val = timeSeries[i] + miter->second[i];

          //Need to be done correctly, and info should be taken from detector geometry
          if (val > numberOfADCCounts) {
            out.str("");
            out << "ADC overflow (ADC = " << val << ", bin " << i << " , station "
                << cIt->GetStationId() << ")\tnoise = " << miter->second[i] << "\tsignal = "
                << timeSeries[i];
            WARNING(out.str());
            val = numberOfADCCounts;
          }
          if (val < 0) {
            out.str("");
            out << "ADC underflow (ADC = " << val << ", bin " << i << " , station "
                << cIt->GetStationId() << ")\tnoise = " << miter->second[i] << "\tsignal = "
                << timeSeries[i];
            WARNING(out.str());
            val = 0;
          }
          timeSeries[i] = val;
        }
      }  // for Channel
    }  //for Station
  }
  return eSuccess;
}

VModule::ResultFlag RdChannelNoiseASCIIImporter::Finish()
{
  // Put any termination or cleanup code here.
  // This method is called once at the end of the run.

  INFO("RdChannelNoiseASCIIImporter::Finish()");

  return eSuccess;
}

string RdChannelNoiseASCIIImporter::GetStringFromStatChan(int statid, int chanid)
{
  ostringstream name;
  name << statid << "_" << chanid;
  return name.str();
}

void RdChannelNoiseASCIIImporter::RemovePedestal(revt::ChannelTimeSeries& timeSeries)
{
  double mean = utl::TraceAlgorithm::Mean(timeSeries, 0, timeSeries.GetSize() - 1);
  int int_mean = int(fabs(mean));
  for (revt::ChannelTimeSeries::Iterator it = timeSeries.Begin(); it != timeSeries.End(); ++it) {
    *it = *it - int_mean;
  }
}

void RdChannelNoiseASCIIImporter::RemovePedestal(std::vector<int>& timeSeries)
{
  double mean = 0.;
  unsigned int N = 0;
  for (std::vector<int>::iterator it = timeSeries.begin(); it != timeSeries.end(); ++it) {
    mean += *it;
    N += 1;
  }
  mean = mean / N;
  int int_mean = int(fabs(mean));
  for (std::vector<int>::iterator it = timeSeries.begin(); it != timeSeries.end(); ++it) {
    *it = *it - int_mean;
  }
}
}