JAdjustZpos.cc

Go to the documentation of this file.

#include <string>
#include <iostream>
#include <iomanip>
#include <vector>
#include <set>
#include <map>
#include <algorithm>
#include <cmath>
#include <sstream>
#include <fstream>
 
#include "TROOT.h"
#include "TFile.h"
 
#include "JAcoustics/JSoundVelocity.hh"
#include "JAcoustics/JEmitter.hh"
#include "JAcoustics/JAcousticsToolkit.hh"
#include "JAcoustics/JHit.hh"
#include "JAcoustics/JFitParameters.hh"
#include "JAcoustics/JGlobalfit.hh"
#include "JAcoustics/JEvent.hh"
#include "JAcoustics/JEvt.hh"
#include "JAcoustics/JEvtToolkit.hh"
#include "JAcoustics/JSupport.hh"
#include "JAcoustics/JTransmission_t.hh"
 
#include "JDetector/JDetector.hh"
#include "JDetector/JDetectorToolkit.hh"
#include "JDetector/JModule.hh"
#include "JDetector/JModuleRouter.hh"
#include "JDetector/JIdealString.hh"
#include "JGeometry3D/JVector3D.hh"
#include "JSupport/JMeta.hh"
 
#include "JSupport/JMultipleFileScanner.hh"
#include "JSupport/JSingleFileScanner.hh"
 
#include "Jeep/JParser.hh"
#include "Jeep/JMessage.hh"
 
int main(int argc, char** argv)
{
  using namespace std;
  using namespace JPP;
 
  string                             inputDetector;
  string                             outputDetector;
  JMultipleFileScanner<JEvent>       inputFile;
  int                                min_working_signals;
  int                                debug;
 
  try {
    JParser<> zap("Adjust z-positions of DOMs without acoustic data and write modified detector");
 
    zap['a'] = make_field(inputDetector,  "input detector file (.datx)");
    zap['o'] = make_field(outputDetector, "output detector file (.datx)");
    zap['f'] = make_field(inputFile,       "output of JAcousticEventBuilder[.sh]");
    zap['m'] = make_field(min_working_signals, "minimum number of acoustic signals to consider a module working") = 10; 
    zap['d'] = make_field(debug,           "debug level") = 0;
    
    zap(argc, argv);
  }
  catch (const exception &e) {
    FATAL(e.what() << endl);
  }
 
  if (inputDetector.empty() || outputDetector.empty() || inputFile.empty()) {
    FATAL("Missing required arguments: -a <input detector> -o <output detector> -f <event files files...>" << endl);
  }
 
  // Load detector
  JDetector detector;
  try {
    load(inputDetector, detector);
  }
  catch (const JException &e) {
    FATAL(e);
  }
 
  JModuleRouter router(detector);
 
  // Build map for acoustic counts per module
  map<int, int> acoustic_counts; // module_id -> count
  for (JDetector::const_iterator m = detector.begin(); m != detector.end(); ++m) {
    acoustic_counts[m->getID()] = 0;
  }
 
  // Scan event files and count acoustic hits per module
  for (const string& file_name : inputFile) {
    try {
      JSingleFileScanner<JEvent> in(file_name, 1);   // read objects one by one
      while (in.hasNext()) {
        const JEvent* evt = in.next();
 
        if (evt->empty()) continue;
 
        // for each transmission/hit in the event, get receiver id and increment count for corresponding module
        for (JEvent::const_iterator hit_it = evt->begin(); hit_it != evt->end(); ++hit_it) {
          const JModuleIdentifier recv_id = hit_it->getID();
          if (router.hasModule(recv_id)) {
            acoustic_counts[recv_id.getID()]++;
          }
        }
      }
    }
    catch (const exception &e) {
      FATAL("Error reading event file " << file_name << ": " << e.what() << endl);
    }
  }
  
  // Map module_id -> has_working based on acoustic counts and cutoff
  map<int, bool> has_working;
  for (JDetector::const_iterator m = detector.begin(); m != detector.end(); ++m) {
    has_working[m->getID()] = (acoustic_counts[m->getID()] >= min_working_signals);
  }
 
  // Group modules by DU
  map<int, vector<const JModule*>> by_du;
  for (JDetector::const_iterator m = detector.begin(); m != detector.end(); ++m) {
    by_du[m->getString()].push_back(&(*m));
  }
 
  // DEBUG: print number of working modules per DU
  if (debug) {
    cout << "Working module counts per DU:" << endl;
    for (const auto &entry : by_du) {
      int du = entry.first;
      const auto &vec = entry.second;
      int working_count = 0;
      for (const auto &mi : vec) {
        if (has_working[mi->getID()]) ++working_count;
      }
      cout << "  DU " << du << ": " << working_count << " working modules out of " << vec.size() << endl;
    }
  }
 
  // Check all DU floor 18 height and classify as ARCA or ORCA based on cutoff
  bool is_arca = false;
  bool is_orca = false;
 
  for (const auto &entry : by_du) {
    int du = entry.first;
    const auto &vec = entry.second;
 
    // Find floor 18 module
    auto it = find_if(vec.begin(), vec.end(),
                      [](const JModule* m){ return m->getFloor() == 18; });
 
    if (it == vec.end()) {
      WARNING("DU " << du << " has no floor 18 module" << endl);
      continue;
    }
 
    double height = (*it)->getPosition().getZ();  
 
    if (height >= STRING_HEIGHT_CUTOFF) {
      is_arca = true;
    } else {
      is_orca = true;
    }
  }
 
  // All DUs must be consistent, otherwise error out
  if (is_arca && is_orca) {
    FATAL("Inconsistent classification: some DUs appear to be ARCA and others ORCA based on floor 18 height" << endl);
  } else if (!is_arca && !is_orca) {
    FATAL("Could not classify as ARCA or ORCA" << endl);
  }
 
  // Select ideal string based on ARCA/ORCA classification
  map<int, double> ideal_string; // floor -> fractional z
  if (is_arca) {
    for (size_t i = 0; i < IDEAL_ARCA_STRING.size(); ++i) {
      ideal_string[static_cast<int>(i)] = IDEAL_ARCA_STRING[i];
    }
  } else {
    for (size_t i = 0; i < IDEAL_ORCA_STRING.size(); ++i) {
      ideal_string[static_cast<int>(i)] = IDEAL_ORCA_STRING[i];
    }
  }
 
  // Debug: print ARCA/ORCA classification and ideal string values
  if (debug) {
    cout << "Detector classified as " << (is_arca ? "ARCA" : "ORCA") << endl;
    cout << "Ideal string fractional z-positions used:" << endl;
    for (const auto &entry : ideal_string) {
      cout << "  Floor " << entry.first << ": " << entry.second << endl;
    }
  }
 
  // Fit per DU: linear model pos_z = frac*L + Z0  => solve normal eqns
  map<int,double> length_fits, offset_fits;
 
  for (const auto &entry : by_du) {
    int du = entry.first;
    const auto &vec = entry.second;
 
    // Filter only working modules that are NOT floor 0
    vector<const JModule*> working_valid;
    for (const JModule* m : vec) {
      if (has_working[m->getID()] && m->getFloor() > 0)
        working_valid.push_back(m);
    }
 
    if (working_valid.size() < 2) continue;
 
    // DEBUG: print all working modules for this DU with their floor and z
    if (debug) {
      cout << "DU " << du << ": Working modules for fit:" << endl;
      for (const auto &mi : working_valid) {
        cout << "  Module ID " << mi->getID() << ", Floor " << mi->getFloor() << ", Z " << FIXED(9,4) << mi->getPosition().getZ() << endl;
      }
    }
 
    double sum_x = 0.0, sum_y = 0.0;
    for (const JModule* m : working_valid) {
      double frac = 0.0;
      auto it = ideal_string.find(m->getFloor());
      if (it != ideal_string.end()) frac = it->second;
      sum_x += frac;
      sum_y += m->getPosition().getZ();
    }
 
    // Calculate means (center data for numerical stability)
    double mean_x = sum_x / working_valid.size();
    double mean_y = sum_y / working_valid.size();
 
    double ss_xx = 0.0, ss_xy = 0.0;
    for (const JModule* m : working_valid) {
      double frac = 0.0;
      auto it = ideal_string.find(m->getFloor());
      if (it != ideal_string.end()) frac = it->second;
        
      // Centering the data
      double diff_x = frac - mean_x;
      ss_xx += diff_x * diff_x;
      ss_xy += diff_x * (m->getPosition().getZ() - mean_y);
    }
 
    // Slope (L) and Intercept (Z0)
    if (std::abs(ss_xx) < 1e-12) continue; // Vertical line protection
 
    double L = ss_xy / ss_xx;
    double Z0 = mean_y - L * mean_x;
 
    length_fits[du] = L;
    offset_fits[du] = Z0;
  }
 
  // DEBUG: print fit results
  if (debug) {
    for (const auto &entry : length_fits) {
      int du = entry.first;
      cout << "DU " << du << ": length fit = " << length_fits[du] << ", offset fit = " << offset_fits[du] << endl;
    }
  }
 
  // Apply adjustments to detector in-memory
  int changed = 0;
  for (JDetector::iterator module = detector.begin(); module != detector.end(); ++module) {
    const int mid = module->getID();
    const int floor = module->getFloor();
    const int du = module->getString();
 
    if (has_working[mid] || floor == 0) continue; // keep original
 
    if (length_fits.count(du) == 0 || offset_fits.count(du) == 0) continue; // no fit
 
    double frac = 0.0;
    auto it = ideal_string.find(floor); 
    if (it != ideal_string.end()) frac = it->second;
 
    double adjusted = frac * length_fits[du] + offset_fits[du];
 
    // keep x,y as before, set new z
    double x = module->getX();
    double y = module->getY();
 
    module->set(JVector3D(x, y, adjusted));
    ++changed;
 
    // print modified module
    cout << "Adjusted module " << mid << " (DU " << du << ", floor " << floor << "): " << adjusted << endl;
  }
 
  // Write updated detector to output file
  detector.comment.add(JMeta(argc, argv));
 
  try {
    store(outputDetector, detector);
  }
  catch (const JException &e) {
    FATAL(e);
  }
 
  cout << "Applied adjustments to " << changed << " modules and wrote " << outputDetector << endl;
 
  return 0;
}