ENH: allow xcorr to work on the component difined in the catalog

apps/scrtdd/descriptions/scrtdd.xml

@@ -164,7 +164,7 @@
                 <description>Maximum data windows lag accepted (secs)..</description>
               </parameter>
               <parameter name="components" type="list:string" default="Z">
-                <description>Priority list of components to be used in cross-correlation. If the component exists for a pair of phases it is used for the cross-correlation, otherwise the next component in the list is tried. Special values are 'H', which computes the L2 norm of the 2 horizontal components and 'T' (and 'R'), which computes the transversal (and radial) component with respect with the event location to which the phase belongs to.</description>
+                <description>Priority list of components to be used in cross-correlation (e.g. Z, N, E, 1, 2, 3). If the component exists for a pair of phases it is used for the cross-correlation, otherwise the next component in the list is tried. Special values are 'H', which computes the L2 norm of the 2 horizontal components and 'T' (and 'R'), which computes the transversal (and radial) component with respect with the event location to which the phase belongs to. If empty the component where the phase was picked is used.</description>
               </parameter>
             </group>
             <group name="s-phase">
@@ -181,7 +181,7 @@
                 <description>Maximum data windows lag accepted (secs).</description>
               </parameter>
               <parameter name="components" type="list:string" default="H">
-                <description>Priority list of components to be used in cross-correlation. If the component exists for a pair of phases it is used for the cross-correlation, otherwise the next component in the list is tried. Special values are 'H', which computes the L2 norm of the 2 horizontal components and 'T' (and 'R'), which computes the transversal (and radial) component with respect with the event location to which the phase belongs to.</description>
+                <description>Priority list of components to be used in cross-correlation (e.g. Z, N, E, 1, 2, 3). If the component exists for a pair of phases it is used for the cross-correlation, otherwise the next component in the list is tried. Special values are 'H', which computes the L2 norm of the 2 horizontal components and 'T' (and 'R'), which computes the transversal (and radial) component with respect with the event location to which the phase belongs to. If empty the component where the phase was picked is used.</description>
               </parameter>
             </group>
             <group name="waveformFiltering">

doc/base/multievent.rst

@@ -155,7 +155,7 @@ E.g. *file phase.csv* ::
 Notes:
 
 * ``type``: mutiple picks are allowed for the same event-station (P,Pn,P1,Pg,S,Sn,S1,Sg), but they must have a different ``type``. However only one P and one S will be used per each event-station (see ``profile.myProfile.catalog.P|S-Phases``).
-* ``channelCode`` is used only for crossCorrelation to select the correct waveform to fetch. However the Orientation Code (the component) of the ``channelCode`` is currently not used (e.g. ``Z`` in ``HHZ``). Instead it is the parameter ``profile.myProfile.crossCorrelation.p|s-phase.components`` that decides the component to use in the crossCorrelation.
+* ``channelCode``: used only for crossCorrelation, it specifies the channel code to use for fetching the waveform. The Orientation Code of the ``channelCode`` (e.g. ``Z`` in ``HHZ``) can be overridden by the parameter ``profile.myProfile.crossCorrelation.p|s-phase.components``.
 * ``lowerUncertainty`` and ``upperUncertainty`` are used only when ``profile.myProfile.solver.aPrioriWeights.usePickUncertainties`` is set to ``true``
 
 With this format it is possible to relocate events that are not stored in any SeisComP database, since all the origins information are contained in those files.
@@ -225,7 +225,8 @@ Relocating a catalog in **"station.csv,event.csv,phase.csv"** file triplet forma
            --inventory-db inventory.xml \
            --verbosity=3 --console=1
 
-The inventory can optionally be empty, which is not an issue if the cross-correlation is not enabled. However if the cross-correlation is used, special waveform transformations (Transversal/Radial component or L2 norm of the horizontal components) become unavailable without an inventory, because those require information contained in the inventory. In this case only existing components can be used.
+The inventory can optionally be empty, which is not an issue if the cross-correlation is not enabled. However when the cross-correlation is used the implication of having an empty invenory is that rtDD doesn't know the orientation of the sensor components and for this reason the special waveform transformations (rotation to the Transversal/Radial component, L2 norm of the horizontal components) become unavailable and should not be selected in ``profile.myProfile.crossCorrelation.p|s-phase.components``. The real componets should be selected instead (e.g. Z, E, N, 1, 2, 3).
+
 This is an **empty inventory**::
 
     <?xml version="1.0" encoding="UTF-8"?>

libs/hdd/dd.cpp

@@ -29,6 +29,7 @@ using Station   = HDD::Catalog::Station;
 using Transform = HDD::Waveform::Processor::Transform;
 using AQ_ACTION = HDD::SolverOptions::AQ_ACTION;
 using HDD::Waveform::getBandAndInstrumentCodes;
+using HDD::Waveform::getOrientationCode;
 
 namespace {
 
@@ -228,6 +229,20 @@ string DD::generateWorkingSubDir(const Event &ev) const
 
 void DD::unloadWaveforms() { createWaveformCache(); }
 
+const std::vector<std::string> DD::xcorrComponents(const Phase &phase) const
+{
+  const auto xcorrCfg = _cfg.xcorr.at(phase.procInfo.type);
+  if (xcorrCfg.components.empty() ||
+      (xcorrCfg.components.size() == 1 && xcorrCfg.components.at(0).empty()))
+  {
+    return {getOrientationCode(phase.channelCode)};
+  }
+  else
+  {
+    return xcorrCfg.components;
+  }
+}
+
 void DD::preloadWaveforms()
 {
   //
@@ -246,11 +261,10 @@ void DD::preloadWaveforms()
         auto eqlrng = _bgCat.getPhases().equal_range(event.id);
         for (auto it = eqlrng.first; it != eqlrng.second; ++it)
         {
-          const Phase &phase  = it->second;
-          TimeWindow tw       = xcorrTimeWindowLong(phase);
-          const auto xcorrCfg = _cfg.xcorr.at(phase.procInfo.type);
+          const Phase &phase = it->second;
+          TimeWindow tw      = xcorrTimeWindowLong(phase);
 
-          for (const string &component : xcorrCfg.components)
+          for (const string &component : xcorrComponents(phase))
           {
             if (func(tw, event, phase, component)) break;
           }
@@ -349,11 +363,10 @@ void DD::dumpWaveforms(const string &basePath)
     auto eqlrng        = _bgCat.getPhases().equal_range(event.id);
     for (auto it = eqlrng.first; it != eqlrng.second; ++it)
     {
-      const Phase &phase  = it->second;
-      TimeWindow tw       = xcorrTimeWindowLong(phase);
-      const auto xcorrCfg = _cfg.xcorr.at(phase.procInfo.type);
+      const Phase &phase = it->second;
+      TimeWindow tw      = xcorrTimeWindowLong(phase);
 
-      for (const string &component : xcorrCfg.components)
+      for (const string &component : xcorrComponents(phase))
       {
         // getAuxProcessor() -> we don't really want to keep the waveforms in
         // memory
@@ -1799,7 +1812,7 @@ DD::preloadNonCatalogWaveforms(Catalog &catalog,
   {
     const Phase &refPhase  = itRef->second;
     const Station &station = catalog.getStations().at(refPhase.stationId);
-    const auto xcorrCfg    = _cfg.xcorr.at(refPhase.procInfo.type);
+    const auto components  = xcorrComponents(refPhase);
 
     // We deal only with real-time event data
     if (refPhase.procInfo.source == Phase::Source::CATALOG) continue;
@@ -1832,7 +1845,7 @@ DD::preloadNonCatalogWaveforms(Catalog &catalog,
         //
         TimeWindow tw = xcorrTimeWindowLong(refPhase);
 
-        for (const string &component : xcorrCfg.components)
+        for (const string &component : components)
         {
           // This doesn't really load the trace but force the request to reach
           // the loader, which will load all the traces later
@@ -2052,14 +2065,12 @@ bool DD::xcorrPhases(const Event &event1,
     }
   }
 
-  auto xcorrCfg  = _cfg.xcorr.at(phase1.procInfo.type);
-  bool performed = false;
-
   //
   // Try all registered components until we are able to perform
   // the cross-correlation
   //
-  for (const string &component : xcorrCfg.components)
+  bool performed = false;
+  for (const string &component : xcorrComponents(phase1))
   {
     performed =
         xcorrPhasesOneComponent(event1, phase1, ph1Cache, event2, phase2,
@@ -2239,7 +2250,11 @@ shared_ptr<const Trace> DD::getWaveform(Waveform::Processor &wfProc,
 
   Catalog::Phase phCopy(ph);
   Transform trans;
-  if (component == "T")
+  if (component.empty())
+  {
+    return nullptr;
+  }
+  else if (component == "T")
   {
     trans = Transform::TRANSVERSAL;
   }

libs/hdd/dd.h

@@ -449,6 +449,9 @@ class DD
 
   void logXCorrSummary(const XCorrCache &xcorr);
 
+  const std::vector<std::string>
+  xcorrComponents(const Catalog::Phase &phase) const;
+
 private:
   const Config _cfg;
 

libs/hdd/waveform.h

@@ -451,7 +451,7 @@ inline std::string getBandAndInstrumentCodes(const std::string &channelCode)
 
 inline std::string getOrientationCode(const std::string &channelCode)
 {
-  if (channelCode.size() == 3) return channelCode.substr(2, 3);
+  if (channelCode.size() == 3) return channelCode.substr(2, 1);
   return "";
 }