/*
Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration
*/
#include "LArRecUtils/LArAutoCorrTotalTool.h"
#include "GaudiKernel/ToolFactory.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/IIncidentSvc.h"
#include "LArElecCalib/LArConditionsException.h"
#include "StoreGate/StoreGateSvc.h"
#include "LArIdentifier/LArOnlineID.h"
#include "LArIdentifier/LArOnline_SuperCellID.h"
#include "LArCabling/LArCablingService.h"
#include "LArCabling/LArSuperCellCablingTool.h"
#include <cmath>
/////////////////////////////////////////////////////////////////////////////
LArAutoCorrTotalTool::LArAutoCorrTotalTool(const std::string& type,
const std::string& name,
const IInterface* parent)
: AthAlgTool(type, name, parent),
m_Nminbias(0),m_MCSym(false),
m_useMixedOFCOpt(false),
m_cablingService(NULL),
m_larmcsym("LArMCSymTool"),m_adc2mevTool("LArADC2MeVTool"),
m_keyShape("LArShape"), m_keyAutoCorr("LArAutoCorr"),
m_keyNoise("LArNoise"), m_keyfSampl("LArfSampl"),
m_keyMinBias("LArMinBias"), m_keyPedestal("LArPedestal"),
m_isMC(true),
m_cacheValid(false),m_loadAtBegin(true),m_deltaBunch(1),m_nsamples(5),m_firstSample(0),m_isSC(false)
{
declareInterface<ILArAutoCorrTotalTool>(this);
declareProperty("keyShape",m_keyShape);
declareProperty("keyAutoCorr",m_keyAutoCorr);
declareProperty("keyNoise",m_keyNoise);
declareProperty("keyPedestal",m_keyPedestal);
declareProperty("keyMinBias",m_keyMinBias);
declareProperty("keyfSampl",m_keyfSampl);
declareProperty("NMinBias",m_Nminbias);
declareProperty("MCSym",m_MCSym);
declareProperty("IsMC",m_isMC);
declareProperty("LoadAtBegin",m_loadAtBegin);
declareProperty("ADC2MeVTool",m_adc2mevTool);
declareProperty("deltaBunch",m_deltaBunch,"Delta between filled bunches in 25 ns units");
declareProperty("NSamples",m_nsamples,"Max number of samples to use");
declareProperty("firstSample",m_firstSample,"First sample to use for in-time event on the full pulse shape");
declareProperty("UseMixedOFCOpt",m_useMixedOFCOpt);
declareProperty("IsSuperCell",m_isSC);
}
/////////////////////////////////////////////////////////////////////////////
StatusCode LArAutoCorrTotalTool::initialize()
{
MsgStream log( msgSvc(), name() );
log << MSG::DEBUG << "LArAutoCorrTotalTool initialize() begin" << endreq;
if ( !m_isSC ) {
const LArOnlineID* laron;
StatusCode sc = detStore()->retrieve(laron,"LArOnlineID");
if (sc.isFailure()) {
log << MSG::ERROR << "Unable to retrieve LArOnlineID from DetectorStore" << endreq;
return StatusCode::FAILURE;
} else m_lar_on_id = (LArOnlineID_Base*) laron;
ToolHandle<LArCablingService> larcab("LArCablingService");
if(larcab.retrieve().isFailure()){
log << MSG::ERROR << "Unable to get CablingService" << endreq;
return StatusCode::FAILURE;
} else m_cablingService = (LArCablingBase*) &(*larcab);
} else {
const LArOnline_SuperCellID* laron;
StatusCode sc = detStore()->retrieve(laron,"LArOnline_SuperCellID");
if (sc.isFailure()) {
log << MSG::ERROR << "Unable to retrieve LArOnlineID from DetectorStore" << endreq;
return StatusCode::FAILURE;
} else m_lar_on_id = (LArOnlineID_Base*) laron;
ToolHandle<LArSuperCellCablingTool> larcab("LArSuperCellCablingTool");
if(larcab.retrieve().isFailure()){
log << MSG::ERROR << "Unable to get CablingService" << endreq;
return StatusCode::FAILURE;
} else m_cablingService = (LArCablingBase*) &(*larcab);
}
//retrieves helpers for LArCalorimeter
m_calo_id_man = CaloIdManager::instance();
if ( m_isSC ) {
m_lar_em_id = m_calo_id_man->getEM_SuperCell_ID();
m_lar_hec_id = m_calo_id_man->getHEC_SuperCell_ID();
m_lar_fcal_id = m_calo_id_man->getFCAL_SuperCell_ID();
} else {
m_lar_em_id = m_calo_id_man->getEM_ID();
m_lar_hec_id = m_calo_id_man->getHEC_ID();
m_lar_fcal_id = m_calo_id_man->getFCAL_ID();
}
if (m_MCSym) {
if (m_larmcsym.retrieve().isFailure()){
log << MSG::ERROR << "Unable to get LArMCSym Tool " << endreq;
return StatusCode::FAILURE;
}
}
if (m_adc2mevTool.retrieve().isFailure())
{
log << MSG::ERROR << "Unable to find tool for LArADC2MeVTool" << endreq;
return StatusCode::FAILURE;
}
else
log << MSG::DEBUG << " -- ILArADC2MeVTool retrieved" << endreq;
if (StatusCode::SUCCESS==detStore()->regFcn(&ILArAutoCorrTotalTool::LoadCalibration,
dynamic_cast<ILArAutoCorrTotalTool*>(this),m_dd_shape,m_keyShape)) {
log << MSG::INFO << "Registered callback for key: " <<m_keyShape << endreq;
} else {
log << MSG::ERROR << "Cannot register testCallback function for key " << m_keyShape << endreq;
}
if (StatusCode::SUCCESS==detStore()->regFcn(&ILArAutoCorrTotalTool::LoadCalibration,
dynamic_cast<ILArAutoCorrTotalTool*>(this),m_dd_autocorr,m_keyAutoCorr)) {
log << MSG::INFO << "Registered callback for key: " << m_keyAutoCorr << endreq;
} else {
log << MSG::ERROR << "Cannot register testCallback function for key " << m_keyAutoCorr << endreq;
}
m_NoPile=false;
if(m_Nminbias<=0) m_NoPile=true;
if(!m_NoPile)
{
if(m_isMC){
if (StatusCode::SUCCESS==detStore()->regFcn(&ILArAutoCorrTotalTool::LoadCalibration,
dynamic_cast<ILArAutoCorrTotalTool*>(this),m_dd_noise,m_keyNoise)) {
log << MSG::INFO << "Registered callback for key: " << m_keyNoise << endreq;
} else {
log << MSG::ERROR << "Cannot register testCallback function for key " << m_keyNoise << endreq;
}
}else{
if (StatusCode::SUCCESS==detStore()->regFcn(&ILArAutoCorrTotalTool::LoadCalibration,
dynamic_cast<ILArAutoCorrTotalTool*>(this),m_dd_pedestal,m_keyPedestal)) {
log << MSG::INFO << "Registered callback for key: " << m_keyPedestal << endreq;
} else {
log << MSG::ERROR << "Cannot register testCallback function for key " << m_keyPedestal << endreq;
}
}
if (StatusCode::SUCCESS==detStore()->regFcn(&ILArAutoCorrTotalTool::LoadCalibration,
dynamic_cast<ILArAutoCorrTotalTool*>(this),m_dd_fSampl,m_keyfSampl)) {
log << MSG::INFO << "Registered callback for key: " << m_keyfSampl << endreq;
} else {
log << MSG::ERROR << "Cannot register testCallback function for key " << m_keyfSampl << endreq;
}
if (StatusCode::SUCCESS==detStore()->regFcn(&ILArAutoCorrTotalTool::LoadCalibration,
dynamic_cast<ILArAutoCorrTotalTool*>(this),m_dd_minbias,m_keyMinBias)) {
log << MSG::INFO << "Registered callback for key: " << m_keyMinBias << endreq;
} else {
log << MSG::ERROR << "Cannot register testCallback function for key " << m_keyMinBias << endreq;
}
// force calling first callback function of LArADC2MeV, and then callback of LArAutoCorrTotalTool
if (StatusCode::SUCCESS==detStore()->regFcn(&ILArADC2MeVTool::LoadCalibration,&(*m_adc2mevTool),
&ILArAutoCorrTotalTool::LoadCalibration,dynamic_cast<ILArAutoCorrTotalTool*>(this))) {
log << MSG::INFO << "Registered callback for LArAutoCorrTool/LArADC2MeVTool" << endreq;
} else {
log << MSG::ERROR << "Cannot register testCallback function for LArAutoCorrTotalTool/LArAdc2MeVTool" << endreq;
}
}
if (m_loadAtBegin) {
log << MSG::DEBUG << "Setting callback function to load calibration at begin of run" << endreq;
// Incident Service:
IIncidentSvc* incSvc;
StatusCode sc = service("IncidentSvc", incSvc);
if (sc.isFailure()) {
log << MSG::ERROR << "Unable to retrieve pointer to IncidentSvc "
<< endreq;
return sc;
}
//start listening to "BeginRun". The incident should be fired AFTER the IOV callbacks and only once.
//const long priority=std::numeric_limits<long>::min(); //Very low priority
// FIXME.. somehow this class can not be a listener.
// incSvc->addListener(this, "BeginRun", priority ,false,true); //single-shot incident
}
//
log << MSG::DEBUG << "LArAutoCorrTotalTool initialize() end" << endreq;
return StatusCode::SUCCESS;
}
/////////////////////////////////////////////////////////////////////////////
StatusCode LArAutoCorrTotalTool::finalize()
{
return StatusCode::SUCCESS;
}
/////////////////////////////////////////////////////////////////////////////
// *** retrieves needed data from the DB ***
StatusCode LArAutoCorrTotalTool::LoadCalibration(IOVSVC_CALLBACK_ARGS_K(keys))
{
MsgStream log( msgSvc(), name() );
log << MSG::DEBUG << "in LoadCalibration " << endreq;
log << MSG::DEBUG << "Callback invoked for " << keys.size() << " keys" << endreq;
m_cacheValid = false;
return StatusCode::SUCCESS;
}
/////////////////////////////////////////////////////////////////////////////
// *** compute some terms of the calculation of the autocorr function ***
StatusCode LArAutoCorrTotalTool::getTerms()
{
MsgStream log( msgSvc(), name() );
log << MSG::DEBUG << "in getAutoCorrTotal" << endreq;
log << MSG::INFO << " Bunch spacing (25 ns units ) " << m_deltaBunch << endreq;
log << MSG::INFO << " N(MB)/bunch crossing " << m_Nminbias << endreq;
// get HWIdentifier iterator
std::vector<HWIdentifier>::const_iterator it =m_lar_on_id->channel_begin();
std::vector<HWIdentifier>::const_iterator it_e =m_lar_on_id->channel_end();
int ngains = (m_isSC ? 1 : 3 );
// resize vector to #(gain) = 3
m_terms.resize(ngains);
int count = 0;
int count2 = 0;
// loop over em Identifiers
log << MSG::DEBUG << "start loop over cells in getAutoCorrTotal" << endreq;
for(;it!=it_e;++it)
{
count ++;
const HWIdentifier id = *it;
unsigned int id32 = id.get_identifier32().get_compact();
if(m_cablingService->isOnlineConnected(id))
{
if(m_MCSym) {
HWIdentifier id2 = m_larmcsym->symOnline(id);
if (id2 != id) continue;
}
count2++;
// mixed OFC optimization: force no pileup in EMB and EMEC-OW
if (m_useMixedOFCOpt) {
const bool isEMB = m_lar_on_id->isEMBchannel(id);
const bool isEMECOW = m_lar_on_id->isEMECOW(id);
if (isEMB || isEMECOW ) {
log << MSG::DEBUG << "No Pileup AutoCorr for ChID 0x" << MSG::hex << id << MSG::dec << endreq;
m_NoPile = true;
} else {
log << MSG::DEBUG << "Using Pileup AutoCorr for ChID 0x" << MSG::hex << id << MSG::dec << endreq;
m_NoPile = false;
}
}
static std::vector<float> empty;
// the Shape is a function of gain
for(int igain=0;igain<ngains;igain++)
{
//:::::::::::::::::::::::::::::::
ILArShape::ShapeRef_t Shape = m_dd_shape->Shape(id,igain);
int nsamples_shape=(int)(Shape.size());
if (nsamples_shape==0)
{
m_terms[igain][id32] = empty;
continue;
}
ILArAutoCorr::AutoCorrRef_t AC = m_dd_autocorr->autoCorr(id,igain);
if (AC.size()==0) {
m_terms[igain][id32] = empty;
continue;
}
m_nsamples_AC_OFC=AC.size()+1;
if (m_nsamples_AC_OFC > m_nsamples) m_nsamples_AC_OFC=m_nsamples;
// fix HEC first sample +1 if the firstSample is 0 and nsamples 4
unsigned int ihecshift=0;
if(m_lar_on_id->isHECchannel(id) && m_nsamples_AC_OFC == 4 && m_firstSample == 0) {
ihecshift=1;
//log << MSG::DEBUG << "Using firstSample +1 for HEC ChID 0x" << MSG::hex << id << MSG::dec << endreq;
}
//:::::::::::::::::::::::::::::::
//NB:
// nsamples_shape = number of samples of the Shape function (e.g 32)
// m_nsamples_AC_OFC = size of AC matrix & OFC vector (e.g 5 in Atlas)
//:::::::::::::::::::::::::::::::
float fSigma2=0.;
if(!m_NoPile)
{
float SigmaNoise;
if(m_isMC)
SigmaNoise = m_dd_noise->noise(id,igain);
else
{
float RMSpedestal =
m_dd_pedestal->pedestalRMS(id,igain);
if(RMSpedestal>(1.0+LArElecCalib::ERRORCODE))
SigmaNoise = RMSpedestal;
else
SigmaNoise = 0.;//(we will have the ERROR message below)
}
float fSampl = m_dd_fSampl->FSAMPL(id);
float MinBiasRMS = m_dd_minbias->minBiasRMS(id);
if(fSampl!=0) MinBiasRMS/=fSampl;
const std::vector<float> *
polynom_adc2mev = &(m_adc2mevTool->ADC2MEV(id,igain));
float Adc2MeV=0.;
if (polynom_adc2mev->size()>0) {
Adc2MeV=(*polynom_adc2mev)[1];
}
if(SigmaNoise!=0 && Adc2MeV!=0)
fSigma2 = pow(MinBiasRMS/(SigmaNoise*Adc2MeV),2);
if(fSampl==0 || SigmaNoise==0 || Adc2MeV==0)
{
if (m_isMC) {
log << MSG::ERROR
<< m_lar_em_id
->show_to_string(m_cablingService->cnvToIdentifier(id))
<< "fSampl ("<<fSampl<<"), SigmaNoise ("
<<SigmaNoise<<") or Adc2MeV ("<<Adc2MeV<<") null "
<<"=> AutoCorrTotal = only AutoCorr elect. part "
<< endreq;
}
fSigma2=0.;
}
//warning: MinBiasRMS is in MeV (at the scale of the hits)
// SigmaNoise is in ADC counts
// so MinBiasRMS/fScale and SigmaNoise*Adc2MeV are the same scale
// (MeV at the scale of the cells)
}//end if m_NoPile
// get in vTerms all the possible non trivial N(N-1)/2 terms of the autocorrelation matrix
int nsize_tot = (m_nsamples_AC_OFC-1)*(m_nsamples_AC_OFC)/2;
std::vector<float> vTerms;
vTerms.resize(2*nsize_tot+m_nsamples_AC_OFC,0.);
//:::::::::::::::::::::::::::::::
for (int j1=0;j1<m_nsamples_AC_OFC-1;j1++) {
for (int j2=j1+1;j2<m_nsamples_AC_OFC;j2++) {
int l=abs(j2-j1)-1;
int index = j1*m_nsamples_AC_OFC- j1*(j1+1)/2 + j2 - (j1+1);
vTerms[index] = AC[l];
}
}
//2nd terms :
for(int j1=0;j1<m_nsamples_AC_OFC-1;++j1) {
for (int j2=j1+1;j2<m_nsamples_AC_OFC;j2++) {
int index = j1*m_nsamples_AC_OFC- j1*(j1+1)/2 + j2 - (j1+1);
float Rij=0;
for(int k=0;k<nsamples_shape;++k) {
if ((j2-j1+k)>=0 && (j2-j1+k)< nsamples_shape) {
int ibunch=0;
if ((j1+m_firstSample+ihecshift-k)%m_deltaBunch == 0) ibunch=1;
Rij += Shape[k] * Shape[j2-j1+k] *ibunch;
}
}
vTerms[nsize_tot+index]= fSigma2*Rij;
}
}
//3rd term : RMS of pileup per samples (multiplied by fSigma2)
for (int j1=0;j1<m_nsamples_AC_OFC;j1++) {
float Rms2i=0;
for(int k=0;k<nsamples_shape;++k) {
int ibunch=0;
if ((j1+m_firstSample+ihecshift-k)%m_deltaBunch == 0) ibunch=1;
Rms2i += pow(Shape[k],2) *ibunch;
}
vTerms[2*nsize_tot+j1]=fSigma2*Rms2i;
}
//storage
m_terms[igain][id32] = vTerms;
}//(loop on gains)
}
else//unconnected
for(unsigned int igain=0;igain<3;igain++)
{
unsigned nsize_tot=(m_nsamples-1)*(m_nsamples)+m_nsamples;
m_terms[igain][id32] = std::vector<float>(nsize_tot, 0.);
}
}
log << MSG::INFO << "LArAutoCorrTotal Ncell " << count << endreq;
log << MSG::INFO << "LArAutoCorrTotal Nsymcell " << count2 << endreq;
log << MSG::DEBUG << "end of loop over cells " << endreq;
m_cacheValid = true;
return StatusCode::SUCCESS;
}
/////////////////////////////////////////////////////////////////////////////
// *** compute AutoCorrTotal (nsamples-1 coeffs) for a given cell ***
const std::vector<double>
LArAutoCorrTotalTool::computeAutoCorr(const std::vector<float>& terms,
float Nminbias) const
{
/*
Rij (elec) + Nminbias*fSigma^2 * sum(k=0->32or800)[G(i-k)G(j-k)]
Rij(total) = ------------------------------------------------------------------
1 + Nminbias*fSigma^2 * sum(k=0->32or800)[G(k)^2]
with: o fSigma = Sigma(E minimum-bias) / Sigma(electronic noise)
o G the Shape function
o Nminbias = number of minimum-bias events per bunch-crossing
(= L/Lo)
if nsamples==5, we need only R01, R02, R03, R04 (for example) (ROO=1)
*/
MsgStream log( msgSvc(), name() );
if(Nminbias<0) Nminbias=m_Nminbias;// takes the value in the property
std::vector<double> vResult;
int tsize = int(sqrt(terms.size()));
int nsize_tot = (tsize-1)*(tsize)/2;
for (int i1=0;i1<tsize-1;i1++) {
for (int i2=i1+1;i2<tsize;i2++) {
int index = i1*tsize - i1*(i1+1)/2 + i2 - (i1+1);
vResult.push_back( (terms[index]+Nminbias*terms[nsize_tot+index])
/ sqrt((1.+Nminbias*terms[2*nsize_tot+i1])*(1.+Nminbias*terms[2*nsize_tot+i2])) );
}
}
return (vResult);
}
//////////////////////////////////////////////////////////////////////////////////////////////////
const std::vector<double>
LArAutoCorrTotalTool::computeRMS(const std::vector<float>& terms,
float Nminbias) const
{
if (Nminbias<0) Nminbias = m_Nminbias;
std::vector<double> vResult;
int tsize = int(sqrt(terms.size()));
vResult.reserve(tsize);
int nsize_tot = (tsize-1)*(tsize)/2;
for (int i=0;i<tsize;i++) {
vResult.push_back(sqrt(1.+Nminbias*terms[2*nsize_tot+i]));
}
return (vResult);
}
/////////////////////////////////////////////////////////////////////////////
// *** retrieve AutoCorrTotal (nsamples*(nsamples-1)/2 coeffs) for a given cell ***
const std::vector<double>
LArAutoCorrTotalTool::autoCorrTotal(const HWIdentifier& CellID,
int gain, float Nminbias) const
{
MsgStream log( msgSvc(), name() );
int thisgain = (m_isSC ? 0 : gain);
if(!m_cacheValid){
LArAutoCorrTotalTool* this2 = const_cast<LArAutoCorrTotalTool*>(this);
StatusCode sc = this2->getTerms();
if (sc.isFailure())
{
log << MSG::ERROR
<< "getTerms failed "
<< endreq;
throw LArConditionsException("Could not compute in LArAutoCorrTotalTool::autoCorrTotal");
}
}
HWIdentifier id;
if (m_MCSym)
id = m_larmcsym->symOnline(CellID);
else
id = CellID;
unsigned int id32 = id.get_identifier32().get_compact();
MAP::const_iterator it = (m_terms[thisgain]).find(id32) ;
if(it == (m_terms[thisgain]).end()){
log << MSG::DEBUG
<< "Unable to find ID = " << CellID << " in m_terms"
<< endreq;
static std::vector<double> empty;
return empty;
}
return ( this->computeAutoCorr((*it).second,Nminbias) );
}
/////////////////////////////////////////////////////////////////////////////
// *** retrieve AutoCorrTotal (4 coeffs) for a given cell ***
const std::vector<double>
LArAutoCorrTotalTool::autoCorrTotal(const Identifier& CellID,
int gain, float Nminbias) const
{
HWIdentifier id = m_cablingService->createSignalChannelID(CellID);
return this->autoCorrTotal(id, gain, Nminbias);
}
void LArAutoCorrTotalTool::handle(const Incident&) {
MsgStream log( msgSvc(), name() );
log << MSG::DEBUG << "In Incident-handle" << endreq;
if(!m_cacheValid){
StatusCode sc = this->getTerms();
if (sc.isFailure())
{
log << MSG::ERROR
<< "getTerms failed "
<< endreq;
throw LArConditionsException("Could not getTerms in LArAutoCorrTotalTool::handle ");
}
}
}
///////////////////////////////////////////////////////////////////////////
const std::vector<double>
LArAutoCorrTotalTool::samplRMS(const HWIdentifier& CellID,
int gain, float Nminbias) const
{
int thisgain = (m_isSC ? 0 : gain);
HWIdentifier id;
if (m_MCSym)
id = m_larmcsym->symOnline(CellID);
else
id = CellID;
unsigned int id32 = id.get_identifier32().get_compact();
MAP::const_iterator it = (m_terms[thisgain]).find(id32) ;
if(it == (m_terms[thisgain]).end()){
MsgStream log( msgSvc(), name() );
log << MSG::ERROR
<< "Unable to find ID = " << CellID << " in m_terms"
<< endreq;
static std::vector<double> empty;
return empty;
}
return ( this->computeRMS((*it).second,Nminbias) );
}
///////////////////////////////////////////////////////////
const std::vector<double>
LArAutoCorrTotalTool::samplRMS(const Identifier& CellID,
int gain, float Nminbias) const
{
HWIdentifier id = m_cablingService->createSignalChannelID(CellID);
return this->samplRMS(id, gain, Nminbias);
}