LuminosityCondAlg.cxx

Go to the documentation of this file.

/*
 * Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration.
 */
#include "LuminosityCondAlg.h"
#include "AthenaPoolUtilities/CondAttrListCollection.h"
#include "StoreGate/ReadCondHandle.h"
#include "StoreGate/WriteCondHandle.h"
#include "CoolKernel/IObject.h"
#include "AthenaKernel/IOVInfiniteRange.h"
#include "CxxUtils/get_unaligned.h"
#include <sstream>
 
 
namespace {
 
// from from BunchCrossingCondAlg, decodes bunch pattern from metadata string
std::vector<float> tokenize( const std::string& pattern ) {
 
  if (pattern == "None") {
    return std::vector<float>(LuminosityCondData::TOTAL_LHC_BCIDS, 1.); 
  }
 
  std::vector< float > result;
  const char* c= pattern.data();
  const char* cEnd= c + pattern.size();
  while(c<cEnd) {
    //This loop swallows actually any string containing numbers
    //separated by non-numbers
    char* end;
    float f=std::strtof(c,&end);
    if (c==end) {//Can't convert, try next
      c++;
    }
    else { //got a value
      result.push_back(f);
      c=end;
    }
  }//end while loop
  return result;
}
 
} // anonymous namespace
 
 
StatusCode
LuminosityCondAlg::initialize()
{
  ATH_CHECK( m_luminosityOutputKey.initialize() );
 
  // May be empty if configured for MC.
  ATH_CHECK( m_onlineLumiCalibrationInputKey.initialize(!m_isMC) );
  ATH_CHECK( m_luminosityFolderInputKey.initialize(!m_isMC) );
 
  // May be empty if configured for data.
  ATH_CHECK( m_mcDigitizationInputKey.initialize(m_isMC) );
  ATH_CHECK( m_eventInfoKey.initialize(m_isMC) );
  ATH_CHECK( m_actualMuKey.initialize(m_isMC) );
  ATH_CHECK( m_averageMuKey.initialize(m_isMC) );
 
  // Only used for run1.
  ATH_CHECK( m_bunchLumisInputKey.initialize(SG::AllowEmpty) );
  ATH_CHECK( m_bunchGroupInputKey.initialize(SG::AllowEmpty) );
  ATH_CHECK( m_fillParamsInputKey.initialize(SG::AllowEmpty) );
  return StatusCode::SUCCESS;
}
 
 
StatusCode
LuminosityCondAlg::execute (const EventContext& ctx) const
{
  auto lumi = std::make_unique<LuminosityCondData>();
  SG::WriteCondHandle<LuminosityCondData> luminosityCondData
    (m_luminosityOutputKey, ctx);
 
  if (m_isMC) {
    // MC case.
    const auto* eventinfo = SG::get(m_eventInfoKey, ctx);
    // define a range for this one LB with one second validity
    EventIDRange range = EventIDRange(EventIDBase(eventinfo->runNumber(),
                                     EventIDBase::UNDEFEVT,
                                     eventinfo->timeStamp(),
                                     eventinfo->timeStampNSOffset(),
                                     eventinfo->lumiBlock()),
                         EventIDBase(eventinfo->runNumber(),
                                     EventIDBase::UNDEFEVT,
                                     eventinfo->timeStamp()+1,
                                     eventinfo->timeStampNSOffset(),
                                     eventinfo->lumiBlock()+1));
 
 
    SG::ReadCondHandle<AthenaAttributeList> digitizationFolder(m_mcDigitizationInputKey, ctx);
    luminosityCondData.addDependency(range);
    
    const float avgMu = eventinfo->averageInteractionsPerCrossing();
    const auto& attr = (**digitizationFolder)[std::string("BeamIntensityPattern")];
    const std::string& sbunches = attr.data<std::string>();
    std::vector<float> bunchpattern = tokenize(sbunches);
    
    if (bunchpattern.size() != LuminosityCondData::TOTAL_LHC_BCIDS) {
      ATH_MSG_ERROR("Decoding MC bunch structure failed, improper number of LHC BCIDs");
      bunchpattern = std::vector<float>(LuminosityCondData::TOTAL_LHC_BCIDS, 1.);
    }
 
    float totintensity = 0.;
    for (size_t i = 0; i < bunchpattern.size(); ++i) {
      totintensity += bunchpattern[i];
      bunchpattern[i] *= (avgMu*m_muToLumi);
    }
 
    lumi->setLbAverageInteractionsPerCrossing(avgMu);
    lumi->setLbAverageLuminosity(totintensity*avgMu*m_muToLumi);
    lumi->setLbLuminosityPerBCIDVector(std::move(bunchpattern));
    lumi->setMuToLumi(m_muToLumi);
  }
  else {
    SG::ReadCondHandle<CondAttrListCollection> luminosityFolder
      (m_luminosityFolderInputKey, ctx);
    luminosityCondData.addDependency(luminosityFolder);
 
    unsigned int preferredChannel;
    unsigned int calibChannel;
    const coral::Blob* bunchInstLumiBlob = nullptr;
    ATH_CHECK( updateAvgLumi (**luminosityFolder,
                              *lumi,
                              preferredChannel,
                              calibChannel,
                              bunchInstLumiBlob) );
 
    ATH_CHECK( updatePerBunchLumi (ctx,
                                   bunchInstLumiBlob,
                                   preferredChannel,
                                   calibChannel,
                                   luminosityCondData,
                                   *lumi) );
  }
 
  
  ATH_CHECK( luminosityCondData.record (std::move (lumi)) );
  return StatusCode::SUCCESS;
}
 
 
StatusCode
LuminosityCondAlg::updateAvgLumi (const CondAttrListCollection& lumiData,
                                  LuminosityCondData& lumi,
                                  unsigned int& preferredChannel,
                                  unsigned int& calibChannel,
                                  const coral::Blob*& bunchInstLumiBlob) const
                                  
{
  preferredChannel = 0;
  calibChannel = 0;
  bunchInstLumiBlob = nullptr;
 
  const coral::AttributeList& attrList = lumiData.attributeList (m_lumiChannel);
  if (attrList.size() == 0 || attrList["Valid"].isNull()) {
    ATH_MSG_DEBUG ("Can't find luminosity information for channel " << m_lumiChannel);
    return StatusCode::SUCCESS;
  }
 
  if (msgLvl (MSG::DEBUG)) {
    std::ostringstream attrStr1;
    attrList.toOutputStream( attrStr1 );
    ATH_MSG_DEBUG( "ChanNum " << m_lumiChannel << " Attribute list "
                   << attrStr1.str() );
  }
 
  // Check data availability
  if (attrList["LBAvInstLumi"].isNull() || attrList["LBAvEvtsPerBX"].isNull()) {
    ATH_MSG_ERROR( " NULL Luminosity information in database " );
    return StatusCode::FAILURE;
  }
 
  // Check validity (don't bother continuing if invalid)
  uint32_t valid = attrList["Valid"].data<cool::UInt32>();
  lumi.setLbAverageValid (valid);
  if (valid & 0x01) {
    if (m_skipInvalid) {
      if (!m_expectInvalid) {
        ATH_MSG_INFO( " Invalid LB Average luminosity ... set lumi to 0" );
      }
      return StatusCode::SUCCESS;
    } else {
      ATH_MSG_DEBUG( " Invalid LB Average luminosity ... continuing because skipInvalid == FALSE" );
    }
  }
 
  // Get preferred channel (needed for per-BCID calculation)
  if (m_lumiChannel == 0u) {
 
    // Check if we have a payload for this (Run2 only)      
    bool hasAlgorithmID = false;
    for (coral::AttributeList::const_iterator attr = attrList.begin();
     attr != attrList.end(); ++attr) {
      if (attr->specification().name() == "AlgorithmID") {
    hasAlgorithmID = true;
    break;
      }
    }
 
    if (hasAlgorithmID) {
      // In Run2, channel 0 should be good.  Leave as is
      preferredChannel = m_lumiChannel;
      calibChannel = attrList["AlgorithmID"].data<cool::UInt32>();
 
    } else {
      // In Run1, we need to recalculate from the actual channel number
      preferredChannel = (valid >> 22);
      calibChannel = preferredChannel;
    }
 
  } else {
    preferredChannel = m_lumiChannel;
    calibChannel = m_lumiChannel;
  }
 
  float LBAvInstLumi = attrList["LBAvInstLumi"].data<cool::Float>();   // Lumi
  float LBAvEvtsPerBX = attrList["LBAvEvtsPerBX"].data<cool::Float>(); // Mu
 
  // Check (and protect for NaN
  if ( std::isnan (LBAvInstLumi) ) {
    ATH_MSG_WARNING( " Luminosity is not a number.. " << LBAvInstLumi << "  ... set it to 0 " );
    LBAvInstLumi=0.;
  }
 
  if ( std::isnan (LBAvEvtsPerBX) ) {
    ATH_MSG_WARNING( " Luminosity is not a number.. " << LBAvEvtsPerBX << "  ... set it to 0 " );
    LBAvEvtsPerBX=0.;
  }
 
  lumi.setLbAverageLuminosity (LBAvInstLumi);
  lumi.setLbAverageInteractionsPerCrossing (LBAvEvtsPerBX);
 
  // Check validity of per-BCID luminosity (will issue warning in recalcPerBCIDLumi
  int perBcidValid = ((valid&0x3ff)/10) % 10;
  if ((perBcidValid > 0) && m_skipInvalid) {
    return StatusCode::SUCCESS;
  }
 
  // Also save per-BCID blob if it exists
  for (coral::AttributeList::const_iterator attr = attrList.begin();
       attr != attrList.end(); ++attr)
  {
    if (attr->specification().name() == "BunchInstLumi") {
      if (!attrList["BunchInstLumi"].isNull())
    bunchInstLumiBlob = &attrList["BunchInstLumi"].data<coral::Blob>();
      break;
    }
  }
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode
LuminosityCondAlg::updatePerBunchLumi (const EventContext& ctx,
                                       const coral::Blob* bunchInstLumiBlob,
                                       unsigned int preferredChannel,
                                       unsigned int calibChannel,
                                       SG::WriteCondHandle<LuminosityCondData>& wHdl,
                                       LuminosityCondData& lumi) const
{
  if (lumi.lbAverageLuminosity() <= 0.) {
    if (!m_expectInvalid) {
      ATH_MSG_INFO( "LBAvInstLumi is zero or negative in updatePerBunchLumi():"
                       << lumi.lbAverageLuminosity());
     wHdl.addDependency(IOVInfiniteRange::infiniteRunLB());
    }
    return StatusCode::SUCCESS;
  }
  bool isValid = true;
  ATH_CHECK( updateMuToLumi (ctx, calibChannel, wHdl, lumi, isValid) );
 
  // Check here if we want to do this the Run1 way (hard) or the Run2 way (easy)
 
  if (!isValid) {
    // Skip if not valid.
  }
  else if (bunchInstLumiBlob != nullptr) { // Run2 way, easy
    ATH_CHECK( updatePerBunchLumiRun2 (*bunchInstLumiBlob,
                                       preferredChannel,
                                       lumi) );
  }
  else { // Run1 way, hard!
    ATH_CHECK( updatePerBunchLumiRun1 (ctx,
                                       preferredChannel,
                                       wHdl,
                                       lumi) );
  }
 
  ATH_MSG_DEBUG( "finished updatePerBunchLumi() for alg: "
                 << preferredChannel );
  return StatusCode::SUCCESS;
}
 
 
StatusCode
LuminosityCondAlg::updateMuToLumi (const EventContext& ctx,
                                   unsigned int calibChannel,
                                   SG::WriteCondHandle<LuminosityCondData>& wHdl,
                                   LuminosityCondData& lumi,
                                   bool& isValid) const
{
  SG::ReadCondHandle<OnlineLumiCalibrationCondData> onlineLumiCalibration
    ( m_onlineLumiCalibrationInputKey, ctx );
  wHdl.addDependency(onlineLumiCalibration);
 
  // This is the only correct way to do this!
  // The division below gives average mu (over all bunches) to total lumi
  float muToLumi = onlineLumiCalibration->getMuToLumi (calibChannel);
 
  // Check if this is reasonable
  if (muToLumi < 0.) {
    ATH_MSG_WARNING(" Found muToLumi = " << muToLumi << " for channel " << calibChannel << ". Try backup channel..." );
    muToLumi = onlineLumiCalibration->getMuToLumi(m_calibBackupChannel);
    ATH_MSG_WARNING(" Found muToLumi = " << muToLumi << " for backup channel " << m_calibBackupChannel);
  }
 
  lumi.setMuToLumi (muToLumi);
 
  // Check validity
  isValid = true;
  int perBcidValid = ((lumi.lbAverageValid()&0x3ff)/10) % 10;
  if ((lumi.lbAverageValid() & 0x03) || (perBcidValid > 0)) {  // Skip if either per-BCID or LBAv is invalid
    isValid = false;
    if (m_skipInvalid) {
      ATH_MSG_WARNING( " Invalid per-BCID luminosity found: "
                       << lumi.lbAverageValid() << "!" );
      return StatusCode::SUCCESS;
    } else {
      ATH_MSG_WARNING( " Invalid per-BCID luminosity found: "
                       << lumi.lbAverageValid()
                       << " continuing because skipInvalid == FALSE" );
    }
  }
 
  // Now check muToLumi and report depending upon whether lumi is valid or not
  if (muToLumi < 0.) {
    if (isValid) {
      ATH_MSG_ERROR(" Found invalid muToLumi = " << muToLumi << " for backup channel " << m_calibBackupChannel << "!");
    } else {
      ATH_MSG_WARNING(" Found invalid muToLumi = " << muToLumi << " for backup channel " << m_calibBackupChannel << "!");
    }
 
    // Don't keep negative values
    muToLumi = 0.;
    lumi.setMuToLumi (muToLumi);
  }
 
  ATH_MSG_DEBUG(" Found muToLumi = " << muToLumi << " for channel "
                << calibChannel );
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode
LuminosityCondAlg::updatePerBunchLumiRun2 (const coral::Blob& bunchInstLumiBlob,
                                           unsigned int preferredChannel,
                                           LuminosityCondData& lumi) const
{
  ATH_MSG_DEBUG( "starting updatePerBunchLumiRun2() for alg: " << preferredChannel );
 
  // Check that the length isn't zero
  if (bunchInstLumiBlob.size() == 0) {
    ATH_MSG_ERROR("BunchInstLumi blob found with zero length!");
    return StatusCode::FAILURE;
  }
 
  // Hardcode the Run2 BLOB decoding (should use CoolLumiUtilities...)
  const uint8_t* ATH_RESTRICT pchar =
    static_cast<const uint8_t*>(bunchInstLumiBlob.startingAddress()); // First byte holds storage size and mode
  unsigned int bss = ((*pchar) % 100) / 10;  // Byte storage size
  unsigned int smod = ((*pchar) % 10);       // Storage mode
  ++pchar;  // Points to next char after header
 
  ATH_MSG_DEBUG( "BunchInstLumi blob found with storage mode " << smod
                 << " and byte storage size " << bss );
 
  // Make sure we have what we think we have
  if (bss != 4 || smod != 1) {
    ATH_MSG_ERROR( "BunchInstLumi blob found with storage mode " << smod << " and byte storage size " << bss << " - Unknown!");
    return StatusCode::FAILURE;
  }
 
  unsigned int nbcids = LuminosityCondData::TOTAL_LHC_BCIDS;
  unsigned int bloblength = bss * nbcids + 1;
 
  if (static_cast<cool::UInt32>(bunchInstLumiBlob.size()) != bloblength) {
    ATH_MSG_ERROR( "BunchRawInstLumi blob found with length"
                   << bunchInstLumiBlob.size() << "in storage mode" << smod
                   <<  ", expecting " << bloblength << "!" );
    return StatusCode::FAILURE;
  }
 
  // Length is correct, read raw data according to packing scheme
  // This is absolute luminosity, so just unpack values into our array
    
  ATH_MSG_DEBUG( "Unpacking lumi value from blob");
  std::vector<float> instLumi (nbcids);
  for (unsigned int i=0; i<nbcids; i++) {
    // Can't use assignment directly because source may be misaligned.
    instLumi[i] = CxxUtils::get_unaligned_float (pchar);
  }
 
  if (msgLvl (MSG::DEBUG)) {
    for (unsigned int i=0; i<nbcids; i++) {
      ATH_MSG_DEBUG( "Bcid: " << i << " Lumi: " << instLumi[i] );
    }
  }
 
  lumi.setLbLuminosityPerBCIDVector (std::move (instLumi));
 
  return StatusCode::SUCCESS;
}
 
 
StatusCode
LuminosityCondAlg::updatePerBunchLumiRun1 (const EventContext& ctx,
                                           unsigned int preferredChannel,
                                           SG::WriteCondHandle<LuminosityCondData>& wHdl,
                                           LuminosityCondData& lumi) const
{
  ATH_MSG_DEBUG( "starting updatePerBunchLumiRun1() for alg: " << preferredChannel );
    
  if (preferredChannel == 0) {
    return StatusCode::SUCCESS;
  }
 
  // Nothing to do if we don't have the ingredients
  if (m_onlineLumiCalibrationInputKey.empty()) {
    ATH_MSG_DEBUG( "OnlineLumiCalibrationInputKey.empty() is TRUE, skipping..." );
    return StatusCode::SUCCESS;
  }
  if (m_bunchLumisInputKey.empty()) {
    ATH_MSG_DEBUG( "BunchLumisInputKey.empty() is TRUE, skipping..." );
    return StatusCode::SUCCESS;
  }
  if (m_bunchGroupInputKey.empty()) {
    ATH_MSG_DEBUG( "BunchGroupTool.empty() is TRUE, skipping..." );
    return StatusCode::SUCCESS;
  }
  if (m_fillParamsInputKey.empty()) {
    ATH_MSG_DEBUG( "FillParamsInputKey.empty() is TRUE, skipping..." );
    return StatusCode::SUCCESS;
  }
 
  SG::ReadCondHandle<OnlineLumiCalibrationCondData> onlineLumiCalibration
    (m_onlineLumiCalibrationInputKey, ctx);
  wHdl.addDependency(onlineLumiCalibration);
  SG::ReadCondHandle<BunchLumisCondData> bunchLumis (m_bunchLumisInputKey, ctx);
  wHdl.addDependency(bunchLumis);
  SG::ReadCondHandle<BunchGroupCondData> bunchGroup (m_bunchGroupInputKey, ctx);
  wHdl.addDependency(bunchGroup);
  SG::ReadCondHandle<FillParamsCondData> fillParams (m_fillParamsInputKey, ctx);
  wHdl.addDependency(fillParams);
 
  const std::vector<unsigned int>& luminousBunches = fillParams->luminousBunches();
  ATH_MSG_DEBUG( "N LuminousBunches:" << luminousBunches.size() );
 
  // Get the raw data for the preferred channel
  const std::vector<float>& rawLumiVec = bunchLumis->rawLuminosity(preferredChannel);
  if (rawLumiVec.empty()) {
    ATH_MSG_DEBUG( "Empty raw luminosity vector" );
    return StatusCode::SUCCESS;
  }
 
  //
  // Calibration step
  //
    
  //  Here we want to go through and calibrate raw values in the luminous bunches only.
  // This is what the OL adds up, and since these are online calibrations, we want to rescale the total
  // to agree to whatever offline tag we are using.
  std::vector<float> calLumiVec (LuminosityCondData::TOTAL_LHC_BCIDS, 0.);
    
  // Update muToLumi while we are at it (also check that calibration exists)
  float muToLumi = onlineLumiCalibration->getMuToLumi (preferredChannel);
  if (muToLumi <= 0.) {
    ATH_MSG_ERROR( " dont have calibration information for preferred channel "
                   << preferredChannel << "!" );
    return StatusCode::FAILURE;
  }
  lumi.setMuToLumi (muToLumi);
 
  double lumiSum = 0.;
  for (unsigned int bcid : luminousBunches) {
    // Don't waste time on zero lumi 
    if (rawLumiVec.at(bcid) <= 0.) {
      ATH_MSG_DEBUG( "Calibrate BCID " << bcid << " with raw "
                     << rawLumiVec.at(bcid) << " -> skipping" );
      continue;
    }
 
    // Calibrate
    if (!onlineLumiCalibration->calibrateLumi(preferredChannel,
                                              rawLumiVec[bcid],
                                              calLumiVec[bcid]))
    {
      ATH_MSG_DEBUG( "Calibrate BCID " << bcid << " with raw " << rawLumiVec[bcid] << " -> calibration failed!" );
      ATH_MSG_WARNING( "Per-BCID calibration failed for bcid " << bcid << " with raw lumi = " << rawLumiVec[bcid] );
      continue;
    }
      
    lumiSum += calLumiVec[bcid];
      
    ATH_MSG_DEBUG( "Calibrate BCID " << bcid << " with raw " << rawLumiVec[bcid] << " -> " << calLumiVec[bcid] );
  }
 
  // Work out scale factor between offline and online estimate
  float offlineOnlineRatio = 1.;
  if (lumiSum > 0.) offlineOnlineRatio = lumi.lbAverageLuminosity() / lumiSum;
    
  ATH_MSG_DEBUG( " Offline/Online scale factor: " << lumi.lbAverageLuminosity()
                 << " / " << lumiSum << " = " << offlineOnlineRatio );
 
  // Make sure we have values for all BCIDs in the physics bunch group
  for (unsigned int bcid : bunchGroup->bunchGroup (1)) {
    // Don't repeat if value already exists
    if (calLumiVec[bcid] > 0.) continue;
    if (rawLumiVec[bcid] <= 0.) continue;
      
    // Calibrate
    if (!onlineLumiCalibration->calibrateLumi(preferredChannel,
                                              rawLumiVec[bcid],
                                              calLumiVec[bcid]))
    {
      ATH_MSG_DEBUG( " -> Calibration failed!" );
      ATH_MSG_WARNING( "Per-BCID calibration failed for bcid " << bcid
                       << " with raw lumi = " << rawLumiVec[bcid] );
      continue;
    }
  }
 
  // Almost done, now we apply the scale factor to all BCIDs
  for (float& lumi : calLumiVec) {
    lumi *= offlineOnlineRatio;
  }
 
  lumi.setLbLuminosityPerBCIDVector (std::move (calLumiVec));
 
  return StatusCode::SUCCESS;
}