CentralTriggerEventBuilder.cc

Go to the documentation of this file.

#include "CentralTriggerEventBuilder.h"

#include <evt/Event.h>

#include <sevt/SEvent.h>
#include <sevt/SEventSimData.h>
#include <sevt/T3.h>
#include <sevt/EventTrigger.h>
#include <sevt/Header.h>
#include <sevt/Station.h>
#include <sevt/StationSimData.h>
#include <sevt/StationTriggerData.h>
#include <sevt/StationGPSData.h>
#include <sevt/PMT.h>
#include <sevt/PMTSimData.h>
#include <sevt/SEventSimData.h>

#include <fevt/FEvent.h>
#include <fevt/Eye.h>
#include <fevt/EyeTriggerData.h>

#include <sdet/SDetector.h>

#include <utl/ErrorLogger.h>
#include <utl/String.h>
#include <utl/UTMPoint.h>
#include <utl/AugerUnits.h>
#include <utl/TabularStream.h>

#include <fwk/CentralConfig.h>

using namespace CentralTriggerEventBuilderNS;
using namespace std;
using namespace fwk;
using namespace utl;


VModule::ResultFlag
CentralTriggerEventBuilder::Init()
{
  fSDetector = &det::Detector::GetInstance().GetSDetector();
  fDefaultWindow = TimeInterval(200*microsecond);
  fDefaultOffset = TimeInterval(0);

  Branch topB = CentralConfig::GetInstance()->GetTopBranch("CentralTriggerEventBuilder");
  topB.GetChild("ForceEventBuildWithoutT3").GetData(fForceEventBuildWithoutT3);

  return eSuccess;
}


VModule::ResultFlag
CentralTriggerEventBuilder::Run(evt::Event& event)
{
  return BuildEvent(event);
}


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


VModule::ResultFlag
CentralTriggerEventBuilder::BuildEvent(evt::Event& event)
{
  if (!event.HasSEvent())
    return eSuccess;
  sevt::SEvent& sEvent = event.GetSEvent();

  // CHECK for FD T3s and count them.
  // REMEMBER the earlierst FD T3 time, as eventual SD trigger time !
  int nTrigEyes = 0;
  TimeStamp t3Time;
  if (event.HasFEvent()) {
    const fevt::FEvent& fEvent = event.GetFEvent();

    for (fevt::FEvent::ConstEyeIterator eyeIt = fEvent.EyesBegin(fevt::ComponentSelector::eHasData);
         eyeIt != fEvent.EyesEnd(fevt::ComponentSelector::eHasData); ++eyeIt)

      if (eyeIt->HasTriggerData() && eyeIt->GetTriggerData().IsT3Accepted()) {
        ++nTrigEyes;
        const TimeStamp& eyeT3time = eyeIt->GetTriggerData().GetT3Time();
        if (!t3Time || eyeT3time < t3Time)
          t3Time = eyeT3time;
      }
  }

  string initAlgo;
  switch (nTrigEyes) {
  case 1:
    initAlgo = "FD";
    break;
  case 2:
    initAlgo = "Stereo";
    break;
  case 3:
    initAlgo = "Trio";
    break;
  case 4:
    initAlgo = "Quad";
    break;
  default:
    initAlgo = "None";
    break;
  }

  // No T3 at all
  if (initAlgo == "None" &&
      !(sEvent.HasSimData() && sEvent.GetSimData().T3sBegin() != sEvent.GetSimData().T3sEnd())) {

    if (!fForceEventBuildWithoutT3) {
      INFO("no FD trigger and no SD T3, continuing");
      return eContinueLoop;
    }

    initAlgo = "Forced";
    INFO("No T3. As specified, forcing build of event in spite of no T3.");
    t3Time = det::Detector::GetInstance().GetTime();

  }

  if (!sEvent.HasTrigger())
    sEvent.MakeTrigger();
  sevt::EventTrigger& evTrig = sEvent.GetTrigger();
  evTrig.SetSender("SSimulation");
  evTrig.SetSDPAngle(0);

  // Decide the time of the event and type/id
  TimeStamp trigTime;
  sevt::T3 t3;

  if (initAlgo != "None") { // Either FD T3 or forced build of event without T3
    evTrig.SetTime(t3Time);
    evTrig.SetAlgorithm(initAlgo);

    // should we really add the entire array for FD T3 ?

    for (sdet::SDetector::StationIterator sIt = fSDetector->StationsBegin();
         sIt != fSDetector->StationsEnd(); ++sIt)
      t3.AddStation(sIt->GetId(), fDefaultOffset, fDefaultWindow);

    ostringstream info;
    info << "Trigger (FD or forced) time set to: " << t3Time;
    INFO(info);
  } else { // SD T3
    const auto& simData = sEvent.GetSimData();
    const int nT3 = simData.GetNumberOfT3s();
    const auto begin = simData.T3sBegin();
    const auto end = simData.T3sEnd();
    {
      ostringstream info;
      info << "Got " << nT3 << " T3 trigger" << (nT3 != 1 ? "s" : "") << ':';
      unsigned int it3 = 0;
      for (auto t3It = begin; t3It != end; ++t3It) {
        info << "\n#" << ++it3 << " n=" << t3It->GetNumberOfStations() << " (" << t3It->GetAlgorithm() << "):\n";
        TabularStream tab("l|l|l");
        tab << "ID" << endc << "offset" << endc << "window" << endr << hline;
        for (auto sIt = t3It->StationsBegin(); sIt != t3It->StationsEnd(); ++sIt)
          tab << sIt->fId << endc << sIt->fOffset/microsecond << endc << sIt->fWindow/microsecond << endr;
        tab << delr;
        info << tab;
      }
      INFO(info);
    }
    auto t3It = begin;
    t3 = *t3It++;
    if (nT3 > 1) {
      INFO("Merging.");
      t3.SetAlgorithm("Merged");
      for ( ; t3It != end; t3It++)
        for (auto sIt = t3.StationsBegin(), end = t3.StationsEnd(); sIt != end; ++sIt)
          t3.AddStation(sIt->fId, sIt->fOffset, sIt->fWindow);
    }
    if (t3) {
      evTrig.SetTime(t3.GetTime());
      evTrig.SetAlgorithm(t3.GetAlgorithm());
      {
        ostringstream info;
        info << "SD Trigger time set to: " << sEvent.GetTrigger().GetTime();
        INFO(info);
      }
    } else
      WARNING("No valid trigger combination found.");
  }

  evTrig.SetId(sEvent.GetHeader().GetId());
  evTrig.SetPreviousId(event.HasFEvent() ? event.GetFEvent().GetHeader().GetId() : 0);

  vector<int> addedTOTStations;
  vector<int> addedMopsStations;
  vector<int> addedRDStations;
  vector<int> addedT1Stations;
  vector<int> addedT2Stations;
  // Iterate over station checking for traces, copying everything to sevt::Station
  int numStations = 0;
  int errorZeroStations = 0;

  bool outOfWindowStations = false;
  TabularStream tab("r|rl");
  tab << "station" << endc << "diff" << endc << " [ns]" << endr << hline;

  for (sevt::T3::StationIterator t3sIt = t3.StationsBegin();
       t3sIt != t3.StationsEnd(); ++t3sIt) {

    const int sId = t3sIt->fId;

    const TimeStamp start = evTrig.GetTime() + (t3sIt->fOffset - t3sIt->fWindow);
    const TimeStamp end = start + 2.*t3sIt->fWindow + TimeInterval(1*microsecond);

    if (!sEvent.HasStation(sId)) {
      ostringstream info;
      info << "skipping station " << sId;
      INFO(info);
      continue;
    }

    sevt::Station& station = sEvent.GetStation(sId);

    if (!station.HasSimData()) {
      ostringstream info;
      info << "station " << station.GetId() << " does not have sim data";
      INFO(info);
      continue;
    }

    // Trigger/GPS data are created and assigned even if they are not T1/T2.
    // If later it happens that there is a T1/T2, the error code is changed
    // The loop continues with the next station after the first T1/T2 found.
    if (!station.HasTriggerData())
      station.MakeTriggerData();
    sevt::StationTriggerData& sTrig = station.GetTriggerData();
    sTrig.SetErrorCode(sevt::StationTriggerData::eT3NotFound);
    sTrig.SetOffsetMicroSecond(t3sIt->fOffset.GetInterval() / microsecond);
    sTrig.SetWindowMicroSecond(t3sIt->fWindow.GetInterval() / microsecond);

    ++numStations;

    // If the station does not have GPS data it shouldn't have triggers,
    // since they are _supposed_ to be created at the same time. However,
    // there is no way to enforce that. It is possible to have a tank trigger
    // simulator without gps simulator, in which case the gps data is created
    // here. Might this cause trouble?
    if (!station.HasGPSData()) {
      //WARNING("Creating GPS data!");
      station.MakeGPSData();
      sevt::StationGPSData& gpsData = station.GetGPSData();

      const int offset = 0;
      const int current100 = 100000000;
      const int next100 = 100000000;
      const int current40 = 40000000;
      const int next40 = 40000000;
      const int previousST = 0;
      const int nextST = 0;
      const int currentST = 0;
      const int second = evTrig.GetTime().GetGPSSecond();
      const int tickFall = 0;
      const int tick =
        (unsigned int)(evTrig.GetTime().GetGPSNanoSecond() /
                       (next100 > 0 ? (1e9 + nextST - currentST)/next100 : 10.0));
      gpsData.SetSecond(second);
      gpsData.SetCurrent100(current100);
      gpsData.SetNext100(next100);
      gpsData.SetCurrent40(current40);
      gpsData.SetNext40(next40);
      gpsData.SetPreviousST(previousST);
      gpsData.SetNextST(nextST);
      gpsData.SetCurrentST(currentST);
      gpsData.SetOffset(offset);
      gpsData.SetTick(tick);
      gpsData.SetTickFall(tickFall);
    }

    sevt::StationSimData& stationSimData = station.GetSimData();

    for (sevt::StationSimData::TriggerTimeIterator tIt = stationSimData.TriggerTimesBegin();
         tIt != stationSimData.TriggerTimesEnd(); ++tIt) {

      const TimeStamp& trigTime = *tIt;

      if (trigTime < start || trigTime > end) {
        outOfWindowStations = true;
        const bool early = (trigTime < start);
        tab << station.GetId() << ' ' << endc
            << ' ' << (early ? (start - trigTime) : (trigTime - end)) / nanosecond << endc
            << (early ? " early" : " late") << endr;
      } else {
        // If it is a T1/T2, copy the traces to the station and set the error code
        const sevt::StationTriggerData& candidateTrigger = stationSimData.GetTriggerData(trigTime);
        if (candidateTrigger.IsT1() || candidateTrigger.IsT2()) {
          ++errorZeroStations;
          sTrig = candidateTrigger;
          sTrig.SetOffsetMicroSecond(t3sIt->fOffset.GetInterval() / microsecond);
          sTrig.SetWindowMicroSecond(t3sIt->fWindow.GetInterval() / microsecond);
          sTrig.SetErrorCode(sevt::StationTriggerData::eNoError);

          station.GetGPSData() = stationSimData.GetGPSData(trigTime);
          for (auto pmtIt = station.PMTsBegin(sdet::PMTConstants::eAnyType);
               pmtIt != station.PMTsEnd(sdet::PMTConstants::eAnyType); ++pmtIt) {

            if (!pmtIt->HasSimData()) {
              ostringstream warn;
              warn << "Station " << sId << ", "
                      "pmt " << pmtIt->GetId() << " does not have sim data but is still in the event?";

              // it will warn for sPMT and scintillator PMT from *accidental stations*
              // as SdAccidentalInjectorKG only uses eWaterCherenkovLarge PMTs.
              // PMTSimData() is made in G4StationSimulatorOG for all PMTs.
              WARNING(warn);
              continue;
            }

            const sevt::PMTSimData& simData = pmtIt->GetSimData();

            // this is needed to avoid NonExistentComponentException in case sPMT and/or scintillator PMT
            // have no FADC traces (no photoelectrons produced)
            if (!simData.HasFADCTrace())
              continue;

            for (const auto gain : { sdet::PMTConstants::eHighGain, sdet::PMTConstants::eLowGain }) {
              for (auto simIt = simData.FADCTriggeredTracesBegin(trigTime, gain),
                     end = simData.FADCTriggeredTracesEnd(trigTime, gain);
                   simIt != end; ++simIt) {
                const auto label = static_cast<sevt::StationConstants::SignalComponent>(simIt->GetLabel());
                if (!pmtIt->HasFADCTrace(label))
                  pmtIt->MakeFADCTrace(label);
                pmtIt->GetFADCTrace(gain, label) = simIt->GetTrace();
              }
              // repeating the whole thing with double-valued traces
              for (auto simIt = simData.FADCTriggeredTracesDBegin(trigTime, gain),
                     end = simData.FADCTriggeredTracesDEnd(trigTime, gain);
                   simIt != end; ++simIt) {
                const auto label = static_cast<sevt::StationConstants::SignalComponent>(simIt->GetLabel());
                if (!pmtIt->HasFADCTrace(label))
                  pmtIt->MakeFADCTrace(label);
                pmtIt->GetFADCTraceD(gain, label) = simIt->GetTrace();
              }
            }

          }

          if (sTrig.IsTimeOverThreshold() || sTrig.IsTimeOverThresholdDeconvoluted())
            addedTOTStations.push_back(sId);
          else if (sTrig.IsMultiplicityOfPositiveSteps())
            addedMopsStations.push_back(sId);
          else if (sTrig.IsRDThreshold())
            addedRDStations.push_back(sId);
          else if (sTrig.IsT2Threshold())
            addedT2Stations.push_back(sId);
          else if (sTrig.IsT1Threshold())
            addedT1Stations.push_back(sId);

          break; // important, this is to take only the first trigger in the time window.
        }
      }
    }
  }

  if (outOfWindowStations) {
    tab << delr;
    cerr << tab << endl;
  }

  // this value includes stations that are part of central trigger
  evTrig.SetNStations(numStations);
  sEvent.SetNErrorZeroStations(errorZeroStations);

  if (!addedTOTStations.empty() ||
      !addedMopsStations.empty() ||
      !addedRDStations.empty() ||
      !addedT1Stations.empty() ||
      !addedT2Stations.empty()) {
    ostringstream info;
    info << "\nAdding stations:";
    if (!addedTOTStations.empty())
      info << '\n' << String::OfSortedIds(addedTOTStations) << " as TOT";
    if (!addedMopsStations.empty())
      info << '\n' << String::OfSortedIds(addedMopsStations) << " as MoPS";
    if (!addedRDStations.empty())
      info << '\n' << String::OfSortedIds(addedRDStations) << " as RD Threshold";
    if (!addedT1Stations.empty())
      info << '\n' << String::OfSortedIds(addedT1Stations) << " as T1 Threshold";
    if (!addedT2Stations.empty())
      info << '\n' << String::OfSortedIds(addedT2Stations) << " as T2 Threshold";
    INFO(info);
  }

  return eSuccess;
}