QwScaler.cc

Go to the documentation of this file.

/**
 * QwScaler.cc
 *
 * Scaler subsystem managing collections of scaler channels with optional
 * normalization, differential operation, and module/channel mapping.
 * Supports SIS3801, STR7200 modules with configurable headers and buffers.
 * Documentation-only edits; runtime behavior unchanged.
 */
 
#include "QwScaler.h"
 
// ROOT headers
#ifdef HAS_RNTUPLE_SUPPORT
#include "ROOT/RNTupleModel.hxx"
#include "ROOT/RField.hxx"
#endif // HAS_RNTUPLE_SUPPORT
 
// Qweak headers
#include "QwParameterFile.h"
 
// Register this subsystem with the factory
REGISTER_SUBSYSTEM_FACTORY(QwScaler);
 
 
/** Define command-line options for scaler subsystem (placeholder). */
void QwScaler::DefineOptions(QwOptions &options)
{
  // Define command line options
}
 
/** Process command-line options for scaler subsystem (placeholder). */
void QwScaler::ProcessOptions(QwOptions &/*options*/)
{
  // Handle command line options
}
 
 
/**
 * Constructor: initialize scaler subsystem with name.
 */
QwScaler::QwScaler(const TString& name)
: VQwSubsystem(name),VQwSubsystemParity(name)
{
  // Nothing, really
}
 
/**
 * Destructor: clean up all owned scaler channel objects.
 */
QwScaler::~QwScaler()
{
  // Delete scalers
  for (size_t i = 0; i < fScaler.size(); i++) {
    delete fScaler.at(i);
  }
  fScaler.clear();
}
 
 
/**
 * Load scaler channel map from file, creating channels and configuring
 * normalization, headers, differential mode, and buffer offsets.
 *
 * @param mapfile Path to the channel map file.
 * @return 0 on success.
 */
Int_t QwScaler::LoadChannelMap(TString mapfile)
{
  // Normalization channel (register default token "1")
  const TString default_norm_channel = "1";
  fName_Map[default_norm_channel] = -1;
  Name_to_Scaler_Map_t::iterator current_norm_channel = fName_Map.find(default_norm_channel);
  double current_norm_factor = 1;
  // Map with normalizations
  std::vector<Name_to_Scaler_Map_t::iterator> norm_channel;
  std::vector<double> norm_factor;
 
  // Include header for this scaler bank
  UInt_t header = 1;
  UInt_t buffer_offset = 0;
 
  // By default the scalers are not differential
  Bool_t differential = false;
 
  // Open the file
  QwParameterFile mapstr(mapfile.Data());
  fDetectorMaps.insert(mapstr.GetParamFileNameContents());
  mapstr.EnableGreediness();
  mapstr.SetCommentChars("!");
  mapstr.AddBreakpointKeyword("norm");
  mapstr.AddBreakpointKeyword("header");
  mapstr.AddBreakpointKeyword("differential");
  mapstr.AddBreakpointKeyword("scaler_buffer_offset");
 
  while (mapstr.ReadNextLine()) {
    RegisterRocBankMarker(mapstr);
 
    mapstr.TrimComment();       // Remove everything after a comment character.
    mapstr.TrimWhitespace();    // Get rid of leading and trailing whitespace
    if (mapstr.LineIsEmpty())  continue;
 
    TString varname, varvalue;
    if (mapstr.HasVariablePair("=", varname, varvalue)) {
      // This is a declaration line.  Decode it.
      varname.ToLower();
      UInt_t value = QwParameterFile::GetUInt(varvalue);
      if (varname == "norm") {
        // Normalization line of format: norm = [ 1 | channel / factor ]
        string dummy = mapstr.GetNextToken("=");
        string channame = mapstr.GetNextToken("/");
        string channorm = mapstr.GetNextToken("/");
        if (fName_Map.count(channame) == 0) {
          // assign a temporarily pointer
          fName_Map[channame] = -1;
        }
        current_norm_channel = fName_Map.find(channame);
        try {
          current_norm_factor = lexical_cast<double>(channorm);
        } catch (boost::bad_lexical_cast&) {
          current_norm_factor = 1;
        }
        QwMessage << "Normalization channel: " << channame << QwLog::endl;
        QwMessage << "Normalization factor: " << current_norm_factor << QwLog::endl;
      } else if (varname == "scaler_buffer_offset") {
        // Buffer offset
        buffer_offset = value;
      } else if (varname == "header") {
        // Header for this block of channels
        header = value;
        QwMessage << "Number of scaler header words: " << header << QwLog::endl;
        if (header > 32)
          QwError << "Is that really what you want?" << QwLog::endl;
      } else if (varname == "differential") {
        // Differential scaler
        try {
          differential = lexical_cast<bool>(varvalue);
        } catch (boost::bad_lexical_cast&) {
          differential = false;
        }
        if (differential)
          QwMessage << "Subsequent scalers will be differential." << QwLog::endl;
        else
          QwMessage << "Subsequent scalers will not be differential." << QwLog::endl;
      }
 
    } else {
 
      //  Break this line into tokens to process it.
      TString modtype = mapstr.GetTypedNextToken<TString>();
      UInt_t  modnum  = mapstr.GetTypedNextToken<UInt_t>();
      UInt_t  channum = mapstr.GetTypedNextToken<UInt_t>();
      TString dettype = mapstr.GetTypedNextToken<TString>();
      TString keyword = mapstr.GetTypedNextToken<TString>();
      modtype.ToUpper();
      dettype.ToUpper();
 
      UInt_t offset = 0;
      if (modtype == "SIS3801" || modtype == "SIS3801D24" || modtype == "SIS3801D32") {
        offset = QwSIS3801D24_Channel::GetBufferOffset(modnum, channum, header)+buffer_offset;
      } else if (modtype == "STR7200") {
        offset = QwSIS3801D32_Channel::GetBufferOffset(modnum, channum, header)+buffer_offset;
      } else {
        QwError << "Unrecognized module type " << modtype << QwLog::endl;
        continue;
      }
 
      // Register data channel type
      Int_t subbank = GetSubbankIndex(fCurrentROC_ID,fCurrentBank_ID);
      if (modnum >= fSubbank_Map[subbank].size())
        fSubbank_Map[subbank].resize(modnum+1);
      if (channum >= fSubbank_Map[subbank].at(modnum).size())
        fSubbank_Map[subbank].at(modnum).resize(channum+1,-1);
      // Add scaler channel
      if (fSubbank_Map[subbank].at(modnum).at(channum) < 0) {
        QwVerbose << "Registering " << modtype << " " << keyword
                  << std::hex
                  << " in ROC 0x" << fCurrentROC_ID << ", bank 0x" << fCurrentBank_ID
                  << std::dec
                  << " at mod " << modnum << ", chan " << channum
                  << QwLog::endl;
        size_t index = fScaler.size();
        VQwScaler_Channel* scaler = 0;
        if (modtype == "SIS3801" || modtype == "SIS3801D24")
          scaler = new QwSIS3801D24_Channel(keyword);
        else if (modtype == "SIS3801D32")
          scaler = new QwSIS3801D32_Channel(keyword);
        else if (modtype == "STR7200")
          scaler = new QwSTR7200_Channel(keyword);
        else {
          QwError << "Unrecognized module type " << modtype << QwLog::endl;
          continue;
        }
 
        // Differential scaler
        scaler->SetDifferentialScaler(differential);
 
        // Register keyword to scaler channel
        fName_Map[keyword] = index;
        fSubbank_Map[subbank].at(modnum).at(channum) = index;
        fModuleChannel_Map[std::pair<Int_t,Int_t>(modnum,channum)] = index;
        // Store scaler
        fScaler.push_back(scaler);
        fNorm.push_back(std::pair<VQwScaler_Channel*,double>(0,1));
        fBufferOffset.push_back(offset);
        // Store current normalization
        norm_channel.push_back(current_norm_channel);
        norm_factor.push_back(current_norm_factor);
      }
    }
  }
 
  // Check for normalization names without actual channels
  for (Name_to_Scaler_Map_t::iterator iter = fName_Map.begin(); iter != fName_Map.end(); iter++) {
    if (iter->second == -1 && iter->first != default_norm_channel) {
      QwError << "Normalization channel " << iter->first << " not found!" << QwLog::endl;
    }
  }
 
  // Fill normalization structure
  fNorm.resize(fScaler.size());
  for (size_t i = 0; i < fScaler.size(); i++) {
    Int_t norm_index = norm_channel.at(i)->second;
    if (norm_index < 0) {
      // No normalization
      fNorm.at(i).first = 0;
      fNorm.at(i).second = 1;
    } else {
      // Normalization
      fNorm.at(i).first = fScaler.at(norm_index);
      fNorm.at(i).second = norm_factor.at(i);
      // Prevent corruption of normalization channel itself
      fNorm.at(norm_index).first = 0;
      fNorm.at(norm_index).second = 1;
    }
  }
 
  mapstr.Close(); // Close the file (ifstream)
  return 0;
}
 
/**
 * Load pedestals and calibration factors for scaler channels.
 *
 * @return 0 on success.
 */
Int_t QwScaler::LoadInputParameters(TString mapfile)
{
  // Open the file
  QwParameterFile mapstr(mapfile.Data());
  fDetectorMaps.insert(mapstr.GetParamFileNameContents());
  while (mapstr.ReadNextLine()) {
    mapstr.TrimComment();    // Remove everything after a comment character
    mapstr.TrimWhitespace(); // Get rid of leading and trailing spaces
    if (mapstr.LineIsEmpty())  continue;
 
    TString varname = mapstr.GetTypedNextToken<TString>();  // name of the channel
    varname.ToLower();
    varname.Remove(TString::kBoth,' ');
    Double_t varped = mapstr.GetTypedNextToken<Double_t>(); // value of the pedestal
    Double_t varcal = mapstr.GetTypedNextToken<Double_t>(); // value of the calibration factor
    QwVerbose << varname << ": " << QwLog::endl;
 
    if (fName_Map.count(varname) != 0) {
      Int_t index = fName_Map[varname];
      QwMessage << "Parameters scaler channel " << varname << ": "
                << "ped = " << varped << ", "
                << "cal = " << varcal << QwLog::endl;
      fScaler[index]->SetPedestal(varped);
      fScaler[index]->SetCalibrationFactor(varcal);
    }
 
  } // end of loop reading all lines of the pedestal file
 
  mapstr.Close(); // Close the file (ifstream)
  return 0;
}
 
/**
 * Clear event data for all scaler channels and reset good event count.
 */
void QwScaler::ClearEventData()
{
  // Clear all scaler channels
  for (size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->ClearEventData();
  }
  // Reset good event count
  fGoodEventCount = 0;
}
 
/**
 * Process the configuration buffer for this subsystem
 * @param roc_id ROC ID
 * @param bank_id Subbank ID
 * @param buffer Buffer to read from
 * @param num_words Number of words left in buffer
 * @return Number of words read
 */
Int_t QwScaler::ProcessConfigurationBuffer(const ROCID_t /*roc_id*/, const BankID_t /*bank_id*/, UInt_t* /*buffer*/, UInt_t /*num_words*/)
{
  return 0;
}
 
/**
 * Process raw buffer data, routing to registered scaler channels.
 *
 * @return Number of words consumed from the buffer.
 */
Int_t QwScaler::ProcessEvBuffer(const ROCID_t roc_id, const BankID_t bank_id, UInt_t* buffer, UInt_t num_words)
{
  // TODO fix :-)
 
  UInt_t words_read = 0;
 
  // Get the subbank index (or -1 when no match)
  Int_t subbank = GetSubbankIndex(roc_id, bank_id);
 
  if (subbank >= 0 && num_words > 0) {
 
    // Read header word
    //UInt_t num_events = buffer[words_read];
    words_read++;
    // TODO Multiscaler functionality
 
    // Read scalers
    for (size_t modnum = 0; modnum < fSubbank_Map[subbank].size(); modnum++) {
      for (size_t channum = 0; channum < fSubbank_Map[subbank].at(modnum).size(); channum++) {
        Int_t index = fSubbank_Map[subbank].at(modnum).at(channum);
        if (index >= 0) {
          UInt_t offset = fBufferOffset.at(index);
          words_read += fScaler[index]->ProcessEvBuffer(&(buffer[offset]), num_words - offset);
        }
      }
    }
    words_read = num_words;
 
  }
 
  return words_read;
}
 
/** Process all scaler channels and apply normalization if configured. */
void QwScaler::ProcessEvent()
{
  // Process the event
  for (size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->ProcessEvent();
  }
 
  // Normalization
  for (size_t i = 0; i < fScaler.size(); i++) {
    if (fNorm.at(i).first) {
      fScaler.at(i)->Scale(fNorm.at(i).second);
      fScaler.at(i)->DivideBy(*(fNorm.at(i).first));
    }
  }
}
 
/** Construct histograms for all scaler channels. */
void QwScaler::ConstructHistograms(TDirectory* folder, TString& prefix)
{
  for(size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->ConstructHistograms(folder, prefix);
  }
}
 
/** Fill histograms for all scaler channels. */
void QwScaler::FillHistograms()
{
  for(size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->FillHistograms();
  }
}
 
/** Construct TTree branches and backing vectors for all scaler channels. */
void QwScaler::ConstructBranchAndVector(TTree *tree, TString & prefix, QwRootTreeBranchVector &values)
{
  for (size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->ConstructBranchAndVector(tree, prefix, values);
  }
}
 
/** Fill tree vector with current scaler channel values. */
void QwScaler::FillTreeVector(QwRootTreeBranchVector &values) const
{
  for(size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->FillTreeVector(values);
  }
}
 
#ifdef HAS_RNTUPLE_SUPPORT
void QwScaler::ConstructNTupleAndVector(std::unique_ptr<ROOT::RNTupleModel>& model, TString& prefix, std::vector<Double_t>& values, std::vector<std::shared_ptr<Double_t>>& fieldPtrs)
{
  for (size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->ConstructNTupleAndVector(model, prefix, values, fieldPtrs);
  }
}
 
void QwScaler::FillNTupleVector(std::vector<Double_t>& values) const
{
  for(size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->FillNTupleVector(values);
  }
}
#endif // HAS_RNTUPLE_SUPPORT
 
/**
 * Assignment operator: copy all scaler channel values.
 */
VQwSubsystem& QwScaler::operator=(VQwSubsystem *value)
{
  if (Compare(value)) {
    VQwSubsystem::operator=(value);
    QwScaler* input = dynamic_cast<QwScaler*>(value);
    for (size_t i = 0; i < fScaler.size(); i++) {
      *(fScaler.at(i)) = *(input->fScaler.at(i));
    }
  }
  return *this; 
}
 
/**
 * Addition-assignment: element-wise addition of scaler channels.
 */
VQwSubsystem& QwScaler::operator+=(VQwSubsystem *value)
{
  if (Compare(value)) {
    QwScaler* input = dynamic_cast<QwScaler*>(value);
    for (size_t i = 0; i < fScaler.size(); i++) {
      *(fScaler.at(i)) += *(input->fScaler.at(i));
    }
  }
  return *this;
}
 
/**
 * Subtraction-assignment: element-wise subtraction of scaler channels.
 */
VQwSubsystem& QwScaler::operator-=(VQwSubsystem *value)
{
  if (Compare(value)) {
    QwScaler* input = dynamic_cast<QwScaler *> (value);
    for (size_t i = 0; i < fScaler.size(); i++) {
      *(fScaler.at(i)) -= *(input->fScaler.at(i));
    }
  }
  return *this;
}
 
 
/**
 * Compute element-wise ratios of scaler channels.
 */
void QwScaler::Ratio(VQwSubsystem *numer, VQwSubsystem *denom)
{
  if (Compare(numer) && Compare(denom)) {
    QwScaler* innumer = dynamic_cast<QwScaler*>(numer);
    QwScaler* indenom = dynamic_cast<QwScaler*>(denom);
    for (size_t i = 0; i < fScaler.size(); i++) {
      fScaler.at(i)->Ratio(*(innumer->fScaler.at(i)), *(indenom->fScaler.at(i)));
    }
  }
}
 
/**
 * Scale all scaler channels by a common factor.
 */
void QwScaler::Scale(Double_t factor)
{
  for (size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->Scale(factor);
  }
}
 
/**
 * Accumulate running sums for all scaler channels.
 */
void QwScaler::AccumulateRunningSum(VQwSubsystem* value, Int_t count, Int_t ErrorMask)
{
  if (Compare(value)) {
    QwScaler* scaler = dynamic_cast<QwScaler*>(value);
    for (size_t i = 0; i < fScaler.size(); i++) {
      fScaler.at(i)->AccumulateRunningSum(scaler->fScaler.at(i), count, ErrorMask);
    }
  }
}
 
/**
 * Deaccumulate running sums for all scaler channels.
 */
void QwScaler::DeaccumulateRunningSum(VQwSubsystem* value, Int_t ErrorMask)
{
  if (Compare(value)) {
    QwScaler* scaler = dynamic_cast<QwScaler*>(value);
    for (size_t i = 0; i < fScaler.size(); i++) {
      fScaler.at(i)->DeaccumulateRunningSum(scaler->fScaler.at(i), ErrorMask);
    }
  }
}
 
/**
 * Calculate running averages for all scaler channels.
 */
void QwScaler::CalculateRunningAverage()
{
  for (size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->CalculateRunningAverage();
  }
}
 
Int_t QwScaler::LoadEventCuts(TString /*filename*/)
{
  return 0;
}
 
Bool_t QwScaler::ApplySingleEventCuts()
{
  return true;
}
 
/*Bool_t QwScaler::CheckForBurpFail(const VQwSubsystem *subsys){
  Bool_t burpstatus = kFALSE;
  VQwSubsystem* tmp = const_cast<VQwSubsystem *>(subsys);
  if(Compare(tmp)) {
    const QwScaler* input = dynamic_cast<const QwScaler*>(subsys);
    for (size_t i = 0; i < input->fScaler.size(); i++) {
      //QwError << "************* test " << typeid(this->fScaler[i]).name() << "*****************" << QwLog::endl;
      burpstatus |= (this->fScaler.at(i))->CheckForBurpFail(input->fScaler.at(i));
    }
  }
  return burpstatus;
}*/
 
void QwScaler::IncrementErrorCounters()
{
}
 
void QwScaler::PrintErrorCounters() const
{
}
 
UInt_t QwScaler::GetEventcutErrorFlag()
{
  return 0;
}
 
Int_t QwScaler::GetChannelIndex(TString channelName, UInt_t module_number)
{
  return 0;
}
 
 
/**
 * Compare two scaler objects
 * @param value Object to compare with
 * @return kTRUE if the object is equal
 */
Bool_t QwScaler::Compare(VQwSubsystem *value)
{
  // Immediately fail on null objects
  if (value == 0) return kFALSE;
 
  // Continue testing on actual object
  Bool_t result = kTRUE;
  if (typeid(*value) != typeid(*this)) {
    result = kFALSE;
 
  } else {
    QwScaler* input = dynamic_cast<QwScaler*> (value);
    if (input->fScaler.size() != fScaler.size()) {
      result = kFALSE;
    }
  }
  return result;
}
 
/**
 * Print some debugging output for the subcomponents
 */
void QwScaler::PrintInfo() const
{
  VQwSubsystemParity::PrintInfo();
 
  QwOut << " there are " << fScaler.size() << " scaler channels" << QwLog::endl;
 
  for (size_t i = 0; i < fScaler.size(); i++) {
    QwOut << " scaler " << i << ": ";
    fScaler.at(i)->PrintInfo();
  }
}
 
/**
 * Print the value for the subcomponents
 */
void QwScaler::PrintValue() const
{
  QwMessage << "=== QwScaler: " << GetName() << " ===" << QwLog::endl;
  for(size_t i = 0; i < fScaler.size(); i++) {
    fScaler.at(i)->PrintValue();
  }
}