/*
 *
 * authors: Ph. Gras (CEA/Saclay), F. Cavallari (INFN/Roma)
 *          some code copied from CalibCalorimetry/EcalTPGTools code
 *          written by P. Paganini and F. Cavallari
 */

#define DB_WRITE_SUPPORT

#include "CalibCalorimetry/EcalSRTools/interface/EcalDccWeightBuilder.h"

#include <limits>
#include <algorithm>
#include <fstream>
#include <iomanip>
#include <TFile.h>
#include <TTree.h>

#include <FWCore/MessageLogger/interface/MessageLogger.h>
#include "SimCalorimetry/EcalSimAlgos/interface/EcalSimParameterMap.h"

#include "DataFormats/EcalDetId/interface/EcalSubdetector.h"

#ifdef DB_WRITE_SUPPORT
#include "OnlineDB/EcalCondDB/interface/EcalCondDBInterface.h"
#include "OnlineDB/EcalCondDB/interface/ODWeightsDat.h"
#endif  //DB_WRITE_SUPPORT defined

#include "CalibCalorimetry/EcalSRTools/src/PasswordReader.h"

using namespace std;
using namespace edm;

const double EcalDccWeightBuilder::weightScale_ = 1024.;

//TODO: handling case of weight encoding saturation: weights shall be downscaled to prevent saturation

EcalDccWeightBuilder::EcalDccWeightBuilder(edm::ParameterSet const& ps)
    : dcc1stSample_(ps.getParameter<int>("dcc1stSample")),
      sampleToSkip_(ps.getParameter<int>("sampleToSkip")),
      nDccWeights_(ps.getParameter<int>("nDccWeights")),
      inputWeights_(ps.getParameter<vector<double> >("inputWeights")),
      mode_(ps.getParameter<string>("mode")),
      dccWeightsWithIntercalib_(ps.getParameter<bool>("dccWeightsWithIntercalib")),
      writeToDB_(ps.getParameter<bool>("writeToDB")),
      writeToAsciiFile_(ps.getParameter<bool>("writeToAsciiFile")),
      writeToRootFile_(ps.getParameter<bool>("writeToRootFile")),
      asciiOutputFileName_(ps.getParameter<string>("asciiOutputFileName")),
      rootOutputFileName_(ps.getParameter<string>("rootOutputFileName")),
      dbSid_(ps.getParameter<string>("dbSid")),
      dbUser_(ps.getParameter<string>("dbUser")),
      dbPassword_(ps.getUntrackedParameter<string>("dbPassword", "")),
      dbTag_(ps.getParameter<string>("dbTag")),
      dbVersion_(ps.getParameter<int>("dbVersion")),
      sqlMode_(ps.getParameter<bool>("sqlMode")),
      geometryToken_(esConsumes()),
      mappingToken_(esConsumes()),
      intercalibConstToken_(esConsumes()),
      calibMap_(emptyCalibMap_),
      ebShape_(consumesCollector()),
      eeShape_(consumesCollector()) {
  if (mode_ == "weightsFromConfig") {
    imode_ = WEIGHTS_FROM_CONFIG;
    if (inputWeights_.size() != (unsigned)nDccWeights_) {
      throw cms::Exception("Config") << "Inconsistent configuration. 'nDccWeights' parameters indicated "
                                     << nDccWeights_ << " weights while parameter 'inputWeights_' contains "
                                     << inputWeights_.size() << " weight values!\n";
    }
  } else if (mode_ == "computeWeights") {
    imode_ = COMPUTE_WEIGHTS;
  } else {
    throw cms::Exception("Config") << "Invalid value ('" << mode_ << "') for parameter mode. "
                                   << "Valid values are: 'weightsFromConfig' and 'computeWeights'\n";
  }
}

void EcalDccWeightBuilder::analyze(const edm::Event& event, const edm::EventSetup& es) {
  const auto mappingHandle = es.getHandle(mappingToken_);
  ecalElectronicsMap_ = mappingHandle.product();

  // Retrieval of intercalib constants
  if (dccWeightsWithIntercalib_) {
    const auto& intercalibConst = es.getData(intercalibConstToken_);
    calibMap_ = intercalibConst.getMap();
  }

  //gets geometry
  geom_ = es.getHandle(geometryToken_);

  //computes the weights:
  computeAllWeights(dccWeightsWithIntercalib_, es);

  //Writing out weights.
  if (writeToAsciiFile_)
    writeWeightToAsciiFile();
  if (writeToRootFile_)
    writeWeightToRootFile();
  if (writeToDB_)
    writeWeightToDB();
}

void EcalDccWeightBuilder::computeAllWeights(bool withIntercalib, const edm::EventSetup& es) {
  const int nw = nDccWeights_;
  int iSkip0_ = sampleToSkip_ >= 0 ? (sampleToSkip_ - dcc1stSample_) : -1;

  EcalSimParameterMap parameterMap;
  const vector<DetId>& ebDetIds = geom_->getValidDetIds(DetId::Ecal, EcalBarrel);
  const vector<DetId>& eeDetIds = geom_->getValidDetIds(DetId::Ecal, EcalEndcap);

  vector<DetId> detIds(ebDetIds.size() + eeDetIds.size());
  copy(ebDetIds.begin(), ebDetIds.end(), detIds.begin());
  copy(eeDetIds.begin(), eeDetIds.end(), detIds.begin() + ebDetIds.size());

  vector<double> baseWeights(nw);  //weight obtained from signal shape
  vector<double> w(nw);            //weight*intercalib
  vector<int> W(nw);               //weight in hw encoding (integrer)
  double prevPhase = numeric_limits<double>::min();

  if (imode_ == WEIGHTS_FROM_CONFIG) {
    assert(inputWeights_.size() == baseWeights.size());
    copy(inputWeights_.begin(), inputWeights_.end(), baseWeights.begin());
  }

  for (vector<DetId>::const_iterator it = detIds.begin(); it != detIds.end(); ++it) {
    double phase = parameterMap.simParameters(*it).timePhase();
    int binOfMax = parameterMap.simParameters(*it).binOfMaximum();

#if 0
    //for debugging...
    edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": ";
    if(it->subdetId()==EcalBarrel){
      edm::LogVerbatim("EcalDccWeightBuilder") << "ieta = " << setw(4) << ((EBDetId)(*it)).ieta()
           << " iphi = " << setw(4) << ((EBDetId)(*it)).iphi() << " ";
    } else if(it->subdetId()==EcalEndcap){
      edm::LogVerbatim("EcalDccWeightBuilder") << "ix = " << setw(3) << ((EEDetId)(*it)).ix()
           << " iy = " << setw(3) << ((EEDetId)(*it)).iy()
           << " iz = " << setw(1) << ((EEDetId)(*it)).iy() << " ";
    } else{
      throw cms::Exception("EcalDccWeightBuilder")
        << "Bug found in " << __FILE__ << ":" << __LINE__ << ": "
        << "Got a detId which is neither tagged as ECAL Barrel "
        << "not ECAL endcap while looping on ECAL cell detIds\n";
    }
    edm::LogVerbatim("EcalDccWeightBuilder") << " -> phase: "  << phase << "\n";
    edm::LogVerbatim("EcalDccWeightBuilder") << " -> binOfMax: " << binOfMax << "\n";
#endif

    try {
      EcalShapeBase* pShape;

      if (it->subdetId() == EcalBarrel) {
        pShape = &ebShape_;
      } else if (it->subdetId() == EcalEndcap) {
        pShape = &eeShape_;
      } else {
        throw cms::Exception("EcalDccWeightBuilder") << "Bug found in " << __FILE__ << ":" << __LINE__ << ": "
                                                     << "Got a detId which is neither tagged as ECAL Barrel "
                                                     << "not ECAL endcap while looping on ECAL cell detIds\n";
      }

      if (phase != prevPhase) {
        if (imode_ == COMPUTE_WEIGHTS) {
          if (it->subdetId() == EcalBarrel) {
            computeWeights(*pShape, binOfMax, phase, dcc1stSample_ - 1, nDccWeights_, iSkip0_, baseWeights);
          }
          prevPhase = phase;
        }
      }
      for (int i = 0; i < nw; ++i) {
        w[i] = baseWeights[i];
        if (withIntercalib)
          w[i] *= intercalib(*it);
      }
      unbiasWeights(w, &W);
      encodedWeights_[*it] = W;
    } catch (std::exception& e) {
      edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": ";
      if (it->subdetId() == EcalBarrel) {
        edm::LogVerbatim("EcalDccWeightBuilder") << "ieta = " << setw(4) << ((EBDetId)(*it)).ieta()
                                                 << " iphi = " << setw(4) << ((EBDetId)(*it)).iphi() << " ";
      } else if (it->subdetId() == EcalEndcap) {
        edm::LogVerbatim("EcalDccWeightBuilder")
            << "ix = " << setw(3) << ((EEDetId)(*it)).ix() << " iy = " << setw(3) << ((EEDetId)(*it)).iy()
            << " iz = " << setw(1) << ((EEDetId)(*it)).iy() << " ";
      } else {
        edm::LogVerbatim("EcalDccWeightBuilder") << "DetId " << (uint32_t)(*it);
      }
      edm::LogVerbatim("EcalDccWeightBuilder") << "phase: " << phase << "\n";
      throw;
    }
  }
}

void EcalDccWeightBuilder::computeWeights(const EcalShapeBase& shape,
                                          int binOfMax,
                                          double timePhase,
                                          int iFirst,
                                          int nWeights,
                                          int iSkip,
                                          vector<double>& result) {
  double sum2 = 0.;
  double sum = 0;
  result.resize(nWeights);

  int nActualWeights = 0;

  const double tzero = -(binOfMax - 1) * 25 + timePhase + shape.timeToRise();  //ns

  for (int i = 0; i < nWeights; ++i) {
    double t_ns = tzero + (iFirst + i) * 25;
    double s = shape(t_ns);
    if (i == iSkip) {
      continue;
    }
    result[i] = s;
    sum += s;
    sum2 += s * s;
    ++nActualWeights;
  }
  for (int i = 0; i < nWeights; ++i) {
    if (i == iSkip) {
      result[i] = 0;
    } else {
      result[i] = (result[i] - sum / nActualWeights) / (sum2 - sum * sum / nActualWeights);
    }
  }
}

int EcalDccWeightBuilder::encodeWeight(double w) { return lround(w * weightScale_); }

double EcalDccWeightBuilder::decodeWeight(int W) { return ((double)W) / weightScale_; }

template <class T>
void EcalDccWeightBuilder::sort(const std::vector<T>& a, std::vector<int>& s, bool decreasingOrder) {
  //performs a bubble sort: adjacent elements are successively swapped 2 by 2
  //until the list is finally sorted.
  bool changed = false;
  s.resize(a.size());
  for (unsigned i = 0; i < a.size(); ++i)
    s[i] = i;
  if (a.empty())
    return;
  do {
    changed = false;
    for (unsigned i = 0; i < a.size() - 1; ++i) {
      const int j = s[i];
      const int nextj = s[i + 1];
      if ((decreasingOrder && (a[j] < a[nextj])) || (!decreasingOrder && (a[j] > a[nextj]))) {
        std::swap(s[i], s[i + 1]);
        changed = true;
      }
    }
  } while (changed);
}

void EcalDccWeightBuilder::unbiasWeights(std::vector<double>& weights, std::vector<int>* encodedWeights) {
  const unsigned nw = weights.size();

  //computes integer weights, weights residuals and weight sum residual:
  vector<double> dw(nw);  //weight residuals due to interger encoding
  vector<int> W(nw);      //integer weights
  int wsum = 0;
  for (unsigned i = 0; i < nw; ++i) {
    W[i] = encodeWeight(weights[i]);
    dw[i] = decodeWeight(W[i]) - weights[i];
    wsum += W[i];
  }

  //sorts weight residuals in decreasing order:
  vector<int> iw(nw);
  sort(dw, iw, true);

  //compensates weight sum residual by adding or substracting 1 to weights
  //starting from:
  // 1) the weight with the minimal signed residual if the correction
  // is positive (wsum<0)
  // 2) the weight with the maximal signed residual if the correction
  // is negative (wsum>0)
  int wsumSign = wsum > 0 ? 1 : -1;
  int i = wsum > 0 ? 0 : (nw - 1);
  while (wsum != 0) {
    W[iw[i]] -= wsumSign;
    wsum -= wsumSign;
    i += wsumSign;
    if (i < 0 || i >= (int)nw) {  //recompute the residuals if a second iteration is
      // needed (in principle, it is not expected with usual input weights), :
      for (unsigned i = 0; i < nw; ++i) {
        dw[i] = decodeWeight(W[i]) - weights[i];
        sort(dw, iw, true);
      }
    }
    if (i < 0)
      i = nw - 1;
    if (i >= (int)nw)
      i = 0;
  }

  //copy result
  if (encodedWeights != nullptr)
    encodedWeights->resize(nw);
  for (unsigned i = 0; i < nw; ++i) {
    weights[i] = decodeWeight(W[i]);
    if (encodedWeights)
      (*encodedWeights)[i] = W[i];
  }
}

double EcalDccWeightBuilder::intercalib(const DetId& detId) {
  // get current intercalibration coeff
  double coef;
  EcalIntercalibConstantMap::const_iterator itCalib = calibMap_.find(detId.rawId());
  if (itCalib != calibMap_.end()) {
    coef = (*itCalib);
  } else {
    coef = 1.;
    edm::LogVerbatim("EcalDccWeightBuilder") << (uint32_t)detId << " not found in EcalIntercalibConstantMap";
  }
#if 0
  edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": ";
  if(detId.subdetId()==EcalBarrel){
    edm::LogVerbatim("EcalDccWeightBuilder") <<  "ieta = " << ((EBDetId)detId).ieta()
	 << " iphi = " << ((EBDetId)detId).iphi();
  } else{
    edm::LogVerbatim("EcalDccWeightBuilder") << "ix = " << ((EEDetId)detId).ix()
	 << " iy = " << ((EEDetId)detId).iy()
	 << " iz = " << ((EEDetId)detId).zside();
  }
  edm::LogVerbatim("EcalDccWeightBuilder") << " coef = " <<  coef << "\n";
#endif
  return coef;
}

void EcalDccWeightBuilder::writeWeightToAsciiFile() {
  string fName = !asciiOutputFileName_.empty() ? asciiOutputFileName_.c_str() : "dccWeights.txt";
  ofstream file(fName.c_str());
  if (!file.good()) {
    throw cms::Exception("Output") << "Failed to open file '" << fName << "'for writing DCC weights\n";
  }

  const char* comment = sqlMode_ ? "-- " : "# ";

  file << comment << "List of weights for amplitude estimation to be used in DCC for\n"
       << comment << "zero suppresssion.\n\n";
  if (!sqlMode_) {
    file << comment << "Note: RU: trigger tower in EB, supercrystal in EE\n"
         << comment << "      xtl: crystal electronic channel id in RU, from 1 to 25\n\n"
         << comment << " DetId    SM  FED RU xtl weights[0..5]...\n";
  }

  if (sqlMode_) {
    file << "variable recid number;\n"
            "exec select COND2CONF_INFO_SQ.NextVal into :recid from DUAL;\n"
            "insert into weights_info (rec_id,tag,version) values (:recid,'"
         << dbTag_ << "'," << dbVersion_ << ");\n";
    file << "\n"
         << comment
         << "index of first sample used in the weighting sum\n"
            "begin\n"
            "  for fedid in "
         << ecalDccFedIdMin << ".." << ecalDccFedIdMax
         << " loop\n"
            "    insert into dcc_weightsample_dat (rec_id, logic_id, sample_id, \n"
            "    weight_number)\n"
            "    values(:recid,fedid,"
         << dcc1stSample_
         << ",1);\n"
            "  end loop;\n"
            "end;\n"
            "/\n";
  } else {
    file << "1st DCC sample: " << dcc1stSample_ << "\n";
  }

  file << "\n" << comment << "list of weights per crystal channel\n";

  for (map<DetId, std::vector<int32_t> >::const_iterator it = encodedWeights_.begin(); it != encodedWeights_.end();
       ++it) {
    const DetId& detId = it->first;

    int fedId;
    int smId;
    int ruId;
    int xtalId;

    //detId ->  fedId, smId, ruId, xtalId
    dbId(detId, fedId, smId, ruId, xtalId);

    char delim = sqlMode_ ? ',' : ' ';

    if (sqlMode_)
      file << "-- detId " << detId.rawId() << "\n"
           << "insert into dcc_weights_dat(rec_id,sm_id,fed_id,"
              "tt_id, cry_id,\n"
              "weight_0,weight_1,weight_2,weight_3,weight_4,weight_5) \n"
              "values ("
              ":recid";

    const vector<int>& weights = it->second;
    if (!sqlMode_)
      file << setw(10) << detId.rawId();
    file << delim << setw(2) << smId;
    file << delim << setw(3) << fedId;
    file << delim << setw(2) << ruId;
    file << delim << setw(2) << xtalId;

    for (unsigned i = 0; i < weights.size(); ++i) {
      file << delim << setw(5) << weights[i];
    }
    if (sqlMode_)
      file << ");";
    file << "\n";
  }
  if (!file.good()) {
    throw cms::Exception("Output") << "Error while writing DCC weights to '" << fName << "' file.";
  }
}
void EcalDccWeightBuilder::writeWeightToRootFile() {
  string fName = !rootOutputFileName_.empty() ? rootOutputFileName_.c_str() : "dccWeights.root";
  TFile file(fName.c_str(), "RECREATE");
  if (file.IsZombie()) {
    throw cms::Exception("Output") << "Failed to open file '" << fName << "'for writing DCC weights\n";
  }
  TTree t("dccWeights", "Weights for DCC ZS filter");
  const int nWeightMax = 20;  //normally n_weights = 6. A different might be used
  //                           used for test purposes.
  struct {
    Int_t detId;
    Int_t fedId;
    Int_t smId;
    Int_t ruId;
    Int_t xtalId;
    Int_t n_weights;
    Int_t weights[nWeightMax];
  } buf;
  t.Branch("weights",
           &buf,
           "rawDetId/I:"
           "feId/I:"
           "smSlotId/I:"
           "ruId/I:"
           "xtalInRuId/I:"
           "n_weights/I:"
           "weights[n_weights]/I");
  for (map<DetId, std::vector<int32_t> >::const_iterator it = encodedWeights_.begin(); it != encodedWeights_.end();
       ++it) {
    buf.detId = it->first.rawId();
    buf.n_weights = it->second.size();

    //detId ->  fedId, smId, ruId, xtalId
    dbId(buf.detId, buf.fedId, buf.smId, buf.ruId, buf.xtalId);

    if (buf.n_weights > nWeightMax) {
      throw cms::Exception("EcalDccWeight") << "Number of weights (" << buf.n_weights << ") for DetId " << buf.detId
                                            << " exceeded maximum limit (" << nWeightMax << ") of root output format. ";
    }
    copy(it->second.begin(), it->second.end(), buf.weights);
    t.Fill();
  }
  t.Write();
  file.Close();
}

#ifndef DB_WRITE_SUPPORT
void EcalDccWeightBuilder::writeWeightToDB() {
  throw cms::Exception("DccWeight") << "Code was compiled without support for writing dcc weights directly "
                                       " into configuration DB. Configurable writeToDB must be set to False. "
                                       "sqlMode can be used to produce an SQL*PLUS script to fill the DB\n";
}
#else   //DB_WRITE_SUPPORT defined
void EcalDccWeightBuilder::writeWeightToDB() {
  edm::LogVerbatim("EcalDccWeightBuilder") << "going to write to the online DB " << dbSid_ << " user " << dbUser_;
  ;
  EcalCondDBInterface* econn;

  try {
    edm::LogVerbatim("EcalDccWeightBuilder") << "Making connection..." << flush;
    const string& filePrefix = string("file:");
    if (dbPassword_.find(filePrefix) == 0) {  //password must be read for a file
      string fileName = dbPassword_.substr(filePrefix.size());
      //substitute dbPassword_ value by the password read from the file
      PasswordReader pr;
      pr.readPassword(fileName, dbUser_, dbPassword_);
    }

    econn = new EcalCondDBInterface(dbSid_, dbUser_, dbPassword_);
    edm::LogVerbatim("EcalDccWeightBuilder") << "Done.";
  } catch (runtime_error& e) {
    edm::LogError("dbconnection") << e.what();
    exit(-1);
  }

  ODFEWeightsInfo weight_info;
  weight_info.setConfigTag(dbTag_);
  weight_info.setVersion(dbVersion_);
  edm::LogVerbatim("EcalDccWeightBuilder") << "Inserting in DB...";

  econn->insertConfigSet(&weight_info);

  int weight_id = weight_info.getId();
  edm::LogVerbatim("EcalDccWeightBuilder") << "WeightInfo inserted with ID " << weight_id;

  vector<ODWeightsDat> datadel;
  datadel.reserve(encodedWeights_.size());

  vector<ODWeightsSamplesDat> dcc1stSampleConfig(nDccs);
  for (int i = ecalDccFedIdMin; i <= ecalDccFedIdMax; ++i) {
    dcc1stSampleConfig[i].setId(weight_id);
    dcc1stSampleConfig[i].setFedId(601 + i);
    dcc1stSampleConfig[i].setSampleId(dcc1stSample_);
    dcc1stSampleConfig[i].setWeightNumber(-1);  //not used.
  }
  econn->insertConfigDataArraySet(dcc1stSampleConfig, &weight_info);

  for (map<DetId, std::vector<int32_t> >::const_iterator it = encodedWeights_.begin(); it != encodedWeights_.end();
       ++it) {
    const DetId& detId = it->first;
    const unsigned nWeights = 6;
    vector<int> weights(nWeights);

    for (unsigned i = 0; i < weights.size(); ++i) {
      //completing the weight vector with zeros in case it has
      //less than 6 elements:
      const vector<int>& w = it->second;
      weights[i] = i < w.size() ? w[i] : 0;
    }

    ODWeightsDat one_dat;
    one_dat.setId(weight_id);

    int fedId;
    int smId;
    int ruId;
    int xtalId;

    //detId ->  fedId, smId, ruId, xtalId
    dbId(detId, fedId, smId, ruId, xtalId);

    one_dat.setSMId(smId);
    one_dat.setFedId(fedId);
    one_dat.setTTId(ruId);
    one_dat.setCrystalId(xtalId);

    one_dat.setWeight0(weights[0]);
    one_dat.setWeight1(weights[1]);
    one_dat.setWeight2(weights[2]);
    one_dat.setWeight3(weights[3]);
    one_dat.setWeight4(weights[4]);
    one_dat.setWeight5(weights[5]);

    datadel.push_back(one_dat);
  }
  econn->insertConfigDataArraySet(datadel, &weight_info);
  edm::LogVerbatim("EcalDccWeightBuilder") << " .. done insertion in DB ";
  delete econn;
  edm::LogVerbatim("EcalDccWeightBuilder") << "closed DB connection ... done";
}
#endif  //DB_WRITE_SUPPORT not defined

void EcalDccWeightBuilder::dbId(const DetId& detId, int& fedId, int& smId, int& ruId, int& xtalId) const {
  const EcalElectronicsId& elecId = ecalElectronicsMap_->getElectronicsId(detId);

  fedId = 600 + elecId.dccId();
  ruId = ecalElectronicsMap_->getElectronicsId(detId).towerId();

  if (detId.subdetId() == EcalBarrel) {
    smId = ((EBDetId)detId).ism();
  } else {
    smId = 10000 - fedId;  //no SM in EE. Use some unique value to satisfy
    //              current DB PK constraints.
  }
  const int stripLength = 5;  //1 strip = 5 crystals in a row
  xtalId = (elecId.stripId() - 1) * stripLength + elecId.xtalId();

#if 0
  edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": FED ID "
       <<  fedId << "\n";

  edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": SM logical ID "
       <<  smId << "\n";
  
  edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": RU ID (TT or SC): "
       <<  ruId << "\n";
  
  edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": strip:"
       <<  elecId.stripId() << "\n";
  
  edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": xtal in strip: "
       <<  elecId.xtalId() << "\n";

  edm::LogVerbatim("EcalDccWeightBuilder") << __FILE__ << ":" << __LINE__ << ": xtalId in RU: "
       <<  xtalId << "\n";
#endif
}