LArBadChannelHunter.cxx File Reference

#include "LArCalibDataQuality/LArBadChannelHunter.h"
#include "LArElecCalib/ILArPedestal.h"
#include "LArRawConditions/LArCaliWaveContainer.h"
#include "LArRawConditions/LArWaveHelper.h"
#include "LArRawConditions/LArWave.h"
#include "LArIdentifier/LArOnlineID.h"
#include "CaloIdentifier/CaloCell_ID.h"
#include "LArRecConditions/LArBadChanBitPacking.h"
#include <fstream>
#include <cmath>

Go to the source code of this file.

Functions
StatusCode LArBadChannelHunter::stop	ATLAS_NOT_THREAD_SAFE ()
	Install fatal handler with default options.

Function Documentation

ATLAS_NOT_THREAD_SAFE()

StatusCode LArBadChannelHunter::stop ATLAS_NOT_THREAD_SAFE ( )

inline

Install fatal handler with default options.

This is meant to be easy to call from python via ctypes.

beginning of the loop of channels

bad bit newly found

known bad bit

for low noisy cells

for high noisy cells

0.01 is used to scale "PER" to the same order of magnitude to "SIG"

smaller deviation: distorted

checking TmaxAmp, Not mixed with MaxAmp and Width

channel information output

Only dead or distorted, or short known BCs are considered below.

index of bc

Definition at line 117 of file LArBadChannelHunter.cxx.

                                                            {
  ATH_CHECK( detStore()->retrieve(m_onlineId, "LArOnlineID") );
  ATH_CHECK( detStore()->retrieve(m_caloId, "CaloCell_ID") );
 
  const ILArPedestal* pedestal = nullptr;
  ATH_CHECK( detStore()->retrieve(pedestal,m_pedKey) );
 
  const LArCaliWaveContainer* caliwave = nullptr;
  ATH_CHECK( detStore()->retrieve(caliwave,m_caliWaveKey) );
 
  if (m_undoCorr) {
    LArCaliWaveContainer* caliwavecomplete= const_cast<LArCaliWaveContainer*>(caliwave);
    ATH_MSG_INFO ( "Undo caliwave corrections now." ) ;
    ATH_CHECK( caliwavecomplete->undoCorrections() );
  }
 
  SG::ReadCondHandle<LArBadChannelCont> readHandle{m_BCKey};
  const LArBadChannelCont *bcCont {*readHandle};
  if(!bcCont) {
     ATH_MSG_ERROR( "Do not have Bad chan container " << m_BCKey.key() );
     return StatusCode::FAILURE;
  }
  SG::ReadCondHandle<LArCalibLineMapping> clHdl{m_CLKey};
  const LArCalibLineMapping *clCont {*clHdl};
  if(!clCont) {
     ATH_MSG_ERROR( "Do not have calib line mapping !!!" );
     return StatusCode::FAILURE;
  }
  SG::ReadCondHandle<LArOnOffIdMapping> cablingHdl{m_cablingKey};
  const LArOnOffIdMapping* cabling=*cablingHdl;
  if(!cabling) {
     ATH_MSG_ERROR( "Do not have cabling object LArOnOffIdMapping" );
     return StatusCode::FAILURE;
  }

  LArWaveHelper waveHelper;
 
  //HWIdentifier lastfeb=HWIdentifier(0);
  //Map of Averages per region (defined by getSymId)
  std::map<unsigned,Average> averageMap;  
  //List of all individual values
  std::vector<CellData> cellData;
 
  std::vector<HWIdentifier>::const_iterator itOnId = m_onlineId->channel_begin();
  std::vector<HWIdentifier>::const_iterator itOnIdEnd = m_onlineId->channel_end();
 
  for(; itOnId!=itOnIdEnd;++itOnId){
    const HWIdentifier chid = *itOnId;
    if (!cabling->isOnlineConnected(chid)) continue;
    //const HWIdentifier febid= m_onlineId->feb_Id(chid);
    float ped=pedestal->pedestalRMS(chid,0);
    const LArCaliWaveVec& cwv=caliwave->get(chid,0);
    if (cwv.size()==0 || ped <=(1.0+LArElecCalib::ERRORCODE)) continue; //Want at least high-gain ped and wave 
    const float ampl=waveHelper.getMaxAmp(cwv[0]);
    const float wid =waveHelper.getWidth(cwv[0])*cwv[0].getDt();  
    const float tmax=waveHelper.getMax(cwv[0])*cwv[0].getDt();
    
    CellData thisCellData;
    thisCellData.m_chid=chid;
    thisCellData.m_ampl=ampl;
    thisCellData.m_wid =wid;
    thisCellData.m_tmax=tmax;
 
    const unsigned regId=getSymId(chid, cabling);
    Average& avreg=averageMap[regId];
 
    for(unsigned igain=CaloGain::LARHIGHGAIN; igain<CaloGain::LARNGAIN; igain++) {
      //get Pedestal RMS
      float ped=pedestal->pedestalRMS(chid,igain);
      if (ped >= (1.0+LArElecCalib::ERRORCODE)) {
    avreg.m_avPedRMS[igain]+=ped;
    avreg.m_avPedRMSSD[igain]+=(ped*ped);
    avreg.m_nPed[igain]++;
    avreg.m_vmedPedRMS[igain].push_back(ped);
    thisCellData.m_pedRMS[igain]=ped;
      }
    }//end loop over gains
 
    avreg.m_avAmpl[0]+=ampl;
    avreg.m_avAmplSD[0]+=ampl*ampl;
    avreg.m_nAmpls[0]++;
    avreg.m_vmedAmpl[0].push_back(ampl);
 
    avreg.m_avWid[0]+=wid;
    avreg.m_avWidSD[0]+=wid*wid;
    avreg.m_nWids[0]++; 
    avreg.m_vmedWid[0].push_back(wid);
 
    avreg.m_avTmax[0]+=tmax;
    avreg.m_avTmaxSD[0]+=tmax*tmax;
    avreg.m_nTmaxs[0]++;
    avreg.m_vmedTmax[0].push_back(tmax);
 
    cellData.push_back(thisCellData);
  }//End loop over cells

  for(auto &[i, average]:averageMap) {
    average.finish(m_recalcPer);
  }
 

  //Keep track to bad cells ...
  typedef LArBadChannel::LArBadChannelEnum::ProblemType PROB_t;
  typedef std::vector<std::pair<HWIdentifier,LArBadChannel> > BCV_t;
  
  BCV_t badChanVec;
  LArBadChanBitPacking packing;
  PROB_t lownoiseProb[CaloGain::LARNGAIN], highnoiseProb[CaloGain::LARNGAIN];
  lownoiseProb[CaloGain::LARHIGHGAIN]=LArBadChannel::LArBadChannelEnum::lowNoiseHGBit; 
  lownoiseProb[CaloGain::LARMEDIUMGAIN]=LArBadChannel::LArBadChannelEnum::lowNoiseMGBit; 
  lownoiseProb[CaloGain::LARLOWGAIN]=LArBadChannel::LArBadChannelEnum::lowNoiseLGBit; 
 
  highnoiseProb[CaloGain::LARHIGHGAIN]=LArBadChannel::LArBadChannelEnum::highNoiseHGBit; 
  highnoiseProb[CaloGain::LARMEDIUMGAIN]=LArBadChannel::LArBadChannelEnum::highNoiseMGBit; 
  highnoiseProb[CaloGain::LARLOWGAIN]=LArBadChannel::LArBadChannelEnum::highNoiseLGBit; 
 
  typedef std::pair<unsigned,std::vector<size_t> > goodAndBad_t;
  typedef std::map<HWIdentifier,goodAndBad_t> goodAndBadMap_t; //the key is the calibLineID
 
  goodAndBadMap_t calibLineMap;
 
  //Second loop: Find bad channels:
  std::vector<CellData>::const_iterator itcells=cellData.begin();
  std::vector<CellData>::const_iterator itcells_e=cellData.end();
 
  for (;itcells!=itcells_e;++itcells) { 
    const CellData& thisCellData=*itcells;
    const HWIdentifier chid = thisCellData.m_chid;
 
    unsigned regId=getSymId(chid, cabling);
    ATH_MSG_VERBOSE ( "Checking " << channelDescription(chid, cabling) ) ;
    std::map<unsigned,Average>::const_iterator febit=averageMap.find(regId);
    if (febit==averageMap.end()) continue;
    const Average& avreg=febit->second;
 
    LArBadChannel problem(0); 
    LArBadChannel bc = bcCont->status(chid); 
 
    for(unsigned igain=CaloGain::LARHIGHGAIN; igain<CaloGain::LARNGAIN; igain++) {
    if (avreg.m_nPed[igain]>2) {//Check if have average ped-rms for this cell & gain
      const float rms=thisCellData.m_pedRMS[igain];
      float lowCut_rms=0.; 
      float higCut_rms=0.; 
      if ( m_cutType=="SIG" ) {
        lowCut_rms=m_lowNoiseTh[igain]*avreg.m_avPedRMSSD[igain];
        higCut_rms=m_highNoiseTh[igain]*avreg.m_avPedRMSSD[igain];
      }  
      else if ( m_cutType=="PER") { 
        lowCut_rms=m_lowNoiseTh[igain]*avreg.m_avPedRMS[igain]*0.01;
        higCut_rms=m_highNoiseTh[igain]*avreg.m_avPedRMS[igain]*0.01;
      }
        
      if (rms==-1)
        ATH_MSG_ERROR ( "No Pedestal found for " << channelDescription(chid,cabling, igain) ) ;
      else {
        ATH_MSG_VERBOSE ( "PedRMS, gain: " << igain << ":" << rms << " Average: " 
                              << avreg.m_avPedRMS[igain] << " Median: " << avreg.m_medPedRMS ) ; 
        if ( (rms-avreg.m_avPedRMS[igain]) > higCut_rms ) {
          packing.setBit(highnoiseProb[igain],problem);
        } 
        else if ( (rms-avreg.m_avPedRMS[igain]) > lowCut_rms ) {
          packing.setBit(lownoiseProb[igain],problem);
        }
        if (fabs(rms-avreg.m_avPedRMS[igain]) > lowCut_rms || (bc.lowNoiseHG()||bc.highNoiseHG()) ) {
          std::string my_status;     
          if (fabs(rms-avreg.m_avPedRMS[igain]) < lowCut_rms && (bc.lowNoiseHG()||bc.highNoiseHG())) 
            my_status="BCN "; 
          else      
            my_status=(bc.lowNoiseHG()||bc.highNoiseHG()) ? "BC " : "NEW ";
 
          ATH_MSG_INFO( my_status << channelDescription(chid,cabling,igain) 
                            << " RMS: " << rms << " ( " << avreg.m_avPedRMS[igain] << " , " 
                            << float(int(10000*(rms-avreg.m_avPedRMS[igain])/avreg.m_avPedRMS[igain]))/100 <<" %) " << ", #Sig: " 
                            << float(int(100*(rms-avreg.m_avPedRMS[igain])/avreg.m_avPedRMSSD[igain]))/100 
                            << " ( " << avreg.m_avPedRMSSD[igain] << " ) " ) ;      
        }// end bad channels information output
      } //end if(rms==-1)
    }//end if have ped rms for this cell & gain
    }//end loop over gains

      if (avreg.m_nAmpls[0]>2 && avreg.m_nWids[0]>2 && avreg.m_nTmaxs[0]>2) {//Check if we have average ampls for this cell & gain
    const float ampl=thisCellData.m_ampl;
    const float wid=thisCellData.m_wid;
    const float tmax=thisCellData.m_tmax;
 
    float lowCut_amp=0.;
    float higCut_amp=0.;
    float lowCut_wid=0.;
    float higCut_wid=0.;
    float Cut_tmax=0.;
 
    if ( m_cutType=="SIG" ) { 
      lowCut_amp=m_distampTh[0]*avreg.m_avAmplSD[0];
      higCut_amp=m_amplTh[0]*avreg.m_avAmplSD[0];
      lowCut_wid=m_distwidTh[0]*avreg.m_avWidSD[0];
      higCut_wid=m_widTh[0]*avreg.m_avWidSD[0];
      Cut_tmax=m_tmaxampTh[0]*avreg.m_avTmaxSD[0];          
    }  
    else if ( m_cutType=="PER") {
      lowCut_amp=m_distampTh[0]*avreg.m_avAmpl[0]*0.01;
      higCut_amp=m_amplTh[0]*avreg.m_avAmpl[0]*0.01;
      lowCut_wid=m_distwidTh[0]*avreg.m_avWid[0]*0.01;
      higCut_wid=m_widTh[0]*avreg.m_avWid[0]*0.01;
      Cut_tmax=m_tmaxampTh[0]*avreg.m_avTmax[0]*0.01;       
    } // end if loop
    
 
    if (ampl==-1 || wid==-1) {
      ATH_MSG_INFO ( "No Amplitude or Width found for " << channelDescription(chid,cabling,0) ) ;
      packing.setBit(LArBadChannel::LArBadChannelEnum::deadReadoutBit,problem);
    }
    else {
      ATH_MSG_VERBOSE ( "Ampl gain: "<< 0<< ":"<< ampl<< " Average: " << avreg.m_avAmpl[0]) ; 
      if (fabs(ampl-avreg.m_avAmpl[0])>higCut_amp || fabs(wid-avreg.m_avWid[0])>higCut_wid) { 
        packing.setBit(LArBadChannel::LArBadChannelEnum::deadReadoutBit,problem);
      } 
      else if (fabs(ampl-avreg.m_avAmpl[0])>lowCut_amp || fabs(wid-avreg.m_avWid[0])>lowCut_wid){ 
            packing.setBit(LArBadChannel::LArBadChannelEnum::distortedBit,problem);
      }
      else if (fabs(tmax-avreg.m_avTmax[0])>Cut_tmax) {
            packing.setBit(LArBadChannel::LArBadChannelEnum::distortedBit,problem);   
//      std::string my_status=(bc.good()) ? "NEW " : "BC ";
//            log << MSG::INFO << my_status << channelDescription(chid,cabling,0) << " Tmax: " << tmax << " ( " << avreg.m_avTmax[0] << " , " 
//          << float(int(10000*(tmax-avreg.m_avTmax[0])/avreg.m_avTmax[0]))/100 << " %) "
//              << " #Sig:" << float(int(100*(tmax-avreg.m_avTmax[0])/avreg.m_avTmaxSD[0]))/100 
//      << " ( " << avreg.m_avTmaxSD[0] << " ) " << endmsg;
      }
//    if (fabs(ampl-avreg.m_avAmpl[0])>lowCut_amp || fabs(wid-avreg.m_avWid[0])>lowCut_wid) {
      if ( problem.deadReadout()||problem.distorted() || bc.deadReadout()||bc.deadCalib()||bc.deadPhys()||bc.distorted()||bc.shortProblem() ){
        std::string my_status;
        if ( !(problem.deadReadout()||problem.distorted())
            &&(bc.deadReadout()||bc.deadCalib()||bc.deadPhys()||bc.distorted()||bc.shortProblem()) )
            my_status="BCN ";
        else    
        my_status=(bc.deadReadout()||bc.deadCalib()||bc.deadPhys()||bc.distorted()||bc.shortProblem()) ? "BC " : "NEW ";
        
        ATH_MSG_INFO ( my_status << channelDescription(chid,cabling,0) 
            << " Amp: " << ampl << " ( " << avreg.m_avAmpl[0] << " , " 
            << float(int(10000*(ampl-avreg.m_avAmpl[0])/avreg.m_avAmpl[0]))/100 << " %) " << " #Sig: " 
        << float(int(100*(ampl-avreg.m_avAmpl[0])/avreg.m_avAmplSD[0]))/100 << " ( " << avreg.m_avAmplSD[0] <<" ) "
            << " FWHM: " << wid << " ( " << avreg.m_avWid[0] << " , " 
        << float(int(10000*(wid-avreg.m_avWid[0])/avreg.m_avWid[0]))/100 << " %) " << " #Sig: " 
        << float(int(100*(wid-avreg.m_avWid[0])/avreg.m_avWidSD[0]))/100 
        << " ( " << avreg.m_avWidSD[0] << " ) "
        << " Tmax: " << tmax << " ( " << avreg.m_avTmax[0] << " , " 
        << float(int(10000*(tmax-avreg.m_avTmax[0])/avreg.m_avTmax[0]))/100 << " %) " << " #Sig:" 
        << float(int(100*(tmax-avreg.m_avTmax[0])/avreg.m_avTmaxSD[0]))/100 
                           ) ;
      }
    }
      }//end if have amplitude for this cell
 
 
      const std::vector<HWIdentifier>& cLids=clCont->calibSlotLine(chid);
      for (const HWIdentifier& hwid : cLids) {
    goodAndBad_t& gb=calibLineMap[hwid];
    if (problem.deadReadout()||problem.distorted())
      gb.second.push_back(badChanVec.size()); 
    else
      ++gb.first;
      }
      if (!problem.good()) badChanVec.push_back(std::make_pair(chid,problem));
 
  }//end loop over channels
 
 
  for(const auto & kv: calibLineMap) {
    const goodAndBad_t& gb=kv.second;
    for (unsigned i=0;i<gb.second.size();i++) {
    if (gb.first==0) {
      ATH_MSG_INFO ( "All channels belonging to calibLine " << channelDescription(badChanVec[gb.second[i]].first, cabling) 
                     << " don't respond to pulses. Assume bad calib line." ) ;
    packing.setBit(LArBadChannel::LArBadChannelEnum::deadReadoutBit,badChanVec[gb.second[i]].second,false);
    packing.setBit(LArBadChannel::LArBadChannelEnum::distortedBit,badChanVec[gb.second[i]].second,false);
    packing.setBit(LArBadChannel::LArBadChannelEnum::deadCalibBit,badChanVec[gb.second[i]].second,true);
      }//end loop over all channels belonging to this calibline
    }
  }//end loop over calibline vector
 
  if (m_outFileName.size()) {
    std::ofstream outfile(m_outFileName.c_str());
    if (!outfile.is_open()) {
      ATH_MSG_ERROR ( "Failed to open output file " << m_outFileName << ". No output will be written." ) ;
    }
    else {
      BCV_t::const_iterator bcvit=badChanVec.begin();
      BCV_t::const_iterator bcvit_e=badChanVec.end();
      for(;bcvit!=bcvit_e;++bcvit) {
    const HWIdentifier chid=bcvit->first;
    const LArBadChannel bc=bcvit->second;
    const HWIdentifier cLid=clCont->calibSlotLine(chid)[0];
    const LArBadChannel bc2=bcCont->status(chid);
    std::string my_ps=(bc2.good())? "NEW " : "BC ";
    if (!bc2.good()&&m_outOnlyNew) continue; // just new
      outfile << m_onlineId->barrel_ec(chid) << " " 
          << m_onlineId->pos_neg(chid) << " "
          << m_onlineId->feedthrough(chid) << " "
          << m_onlineId->slot(chid) << " "
          << m_onlineId->channel(chid) << " "
          << m_onlineId->channel(cLid) << " "
          << packing.stringStatus(bc) << " "
          << my_ps << std::endl;
//        << packing.stringStatus(bc) << std::endl;
      }
      outfile.close();
    } //end if out file open
  } // end if have out file name
  return StatusCode::SUCCESS; 
}