QChannelCovariance.cc

Go to the documentation of this file.

#include "QChannelCovariance.hh"
#include <sstream>
#include <cmath>
 
QObjectImp(QChannelCovariance);
 
bool CompareCorrelations(const std::pair<int,double>& lhs, const std::pair<int,double>& rhs)
{
    return lhs.second > rhs.second;
}
 
 
QChannelCovariance::QChannelCovariance()
{
}
 
QChannelCovariance::QChannelCovariance(size_t nfreq, const std::vector<int> &channels)
{
  SetNFreq(nfreq);
  SetChannels(channels);
}
 
void QChannelCovariance::SetNFreq(size_t nfreq)
{
  Clear();
  fCovariance.resize(nfreq);
}
 
void QChannelCovariance::SetChannels(const std::vector<int> &channels)
{
  for(size_t f = 0; f < fCovariance.size(); f++) {
      fCovariance[f].Resize(channels.size(),channels.size());
      fCovariance[f].Initialize(0,0);
  }
  for(size_t i = 0; i < channels.size(); i++) {
      fChannels[channels[i]]=i;
  }
  fNevents.Resize(channels.size(),channels.size());
  fNevents.Initialize(0);
}
 
void QChannelCovariance::SetCovariance(int chI, int chJ,const Diana::QVectorC& cov) 
{
    if(cov.Size() != fCovariance.size()) {
        QError err = QERR_SIZE_NOT_MATCH;
        std::stringstream msg;
        msg<<"Size of given QVectorC, "
        <<cov.Size()<<", is wrong; expects" <<fCovariance.size();
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
        return;
    }
 
    int i,j;
    Channels2Indices(chI,chJ,i,j);
 
    for(size_t f = 0; f < fCovariance.size(); f++) {
        fCovariance[f](i,j) = cov[f];
    }
}
 
void QChannelCovariance::SetNEvents(int chI, int chJ, int nev)
{
    int i,j;
    Channels2Indices(chI,chJ,i,j);
    fNevents(i,j) = nev; 
}
 
std::vector<int> QChannelCovariance::GetChannels() const
{
  std::map<int,int>::const_iterator iter = fChannels.begin();
  std::vector<int> res(fChannels.size());
  for(size_t i = 0; iter != fChannels.end(); i++) {
      res[i]= iter->first;
      iter++;
  }
  return res;
}
 
 
Diana::QVectorC QChannelCovariance::GetCovariance(int chI, int chJ) const
{
    int i,j;
    Channels2Indices(chI,chJ,i,j);
 
    Diana::QVectorC res(fCovariance.size());
    for(size_t f = 0; f < fCovariance.size(); f++) {
        res[f] = fCovariance[f](i,j);
    }
    return res;
}
 
Diana::QVectorC QChannelCovariance::GetCovarianceByIndex(int i, int j) const
{
    Diana::QVectorC res(fCovariance.size());
    for(size_t f = 0; f < fCovariance.size(); f++) {
        res[f] = fCovariance[f](i,j);
    }
    return res;
}
 
 
 
int QChannelCovariance::GetNEvents(int chI, int chJ) const
{
    int i,j;
    Channels2Indices(chI,chJ,i,j);
    return (int)fNevents(i,j);
}
 
QChannelCovariance QChannelCovariance::Inverse() const
{
    QChannelCovariance res(*this);
    for(size_t f = 1; f < fCovariance.size(); f++) {
        res.fCovariance[f].Invert();
    }
    return res;
}
 
QChannelCovariance QChannelCovariance::GetSubMatrix(const std::vector<int> &channels) const
{
    QChannelCovariance res(fCovariance.size(),channels);
    for(size_t i = 0; i < channels.size(); i++) {
        for(size_t j = 0; j < channels.size(); j++) {
            const int chanI = channels[i];
            const int chanJ = channels[j];
            res.SetCovariance(chanI,chanJ,GetCovariance(chanI,chanJ));
            res.SetNEvents(chanI,chanJ,GetNEvents(chanI,chanJ));
        }
    }
    return res;
}
 
 
 
 
Diana::QVectorC QChannelCovariance::GetCorrelation(int chI, int chJ) const
{
    Diana::QVectorC covII, covJJ, covIJ;
    covII = GetCovariance(chI,chI);
    covJJ = GetCovariance(chJ,chJ);
    covIJ = GetCovariance(chI,chJ);
 
    Diana::QVectorC res(covII.Size());
    res[0].SetRe(0);
    res[0].SetIm(0);
    for(size_t f = 1; f < covII.Size(); f++) {
        const double div = sqrt(fabs(covII[f].Re()*covJJ[f].Re()));
        res[f] = covIJ[f]/div;
    }
    return res;
}
 
double QChannelCovariance::GetMeanCorrelation(int chI, int chJ, int fmin, int fmax) const
{
    Diana::QVectorC covII = GetCovariance(chI,chI);
    Diana::QVector rho2 = GetCorrelation(chI,chJ).GetModulusSquare();
 
    double num = 0., den = 0.;
    if(fmin <= 0 || fmin >= (int)rho2.Size()) fmin = 1;
    if(fmax < 0 || fmax >= (int)rho2.Size()) fmax = -1+rho2.Size();
    
    for(int f = fmin; f < fmax; f++) {
        num += rho2[f]*fabs(covII[f].Re());
        den += fabs(covII[f].Re());
    }
    double res = 0.;
    if(den != 0) res = sqrt(num/den);
    return res;
}
 
std::vector<int> QChannelCovariance::GetChannelsByCorrelation(int chI, int fmin, int fmax) const
{
    std::vector<std::pair<int,double> > pairs;
    std::vector<int> ch = GetChannels();
    for(size_t i  = 0; i < ch.size();  i++) {
        int chJ = ch[i];
        if(chJ != chI) {
            double corr = GetMeanCorrelation(chI,chJ,fmin,fmax);
            std::pair<int,double> p(chJ,corr);
            pairs.push_back(p);
        }
    }
    std::sort(pairs.begin(),pairs.end(),CompareCorrelations);
    std::vector<int> res; 
    for(size_t i  = 0; i < pairs.size();  i++) {
        res.push_back(pairs[i].first);
    }
    return res;
}
 
 
 
void QChannelCovariance::Dump(std::ostream& o) const
{
    std::vector<int> ch = GetChannels();
    o<<"# of chans: "<<ch.size()<<std::endl;
    o<<"# of freqs: "<<fCovariance.size();
 
    if(ch.empty()) return;
    o<<std::endl;
    if(ch.size() == 1) {
        o<<"channel "<<ch[0];
        if(fCovariance.size()<2) return;
        o<<": Covariance[f=1] = "<<fCovariance[1](ch[0],ch[0]);
        o<<", Nevents = "<<fNevents(ch[0],ch[0]);
        return;
    } else {
        o<<"first 2 chans ("<<ch[0]<<","<<ch[1]<<")";
        if(fCovariance.size()<2) return;
        o<<": Covariance[f=1] = "<<fCovariance[1](ch[0],ch[1]);
        o<<", Nevents = "<<fNevents(ch[0],ch[0]);
    }
    return;
 
}
 
void QChannelCovariance::DumpMeanCorrelationMatrix(std::ostream& o,const int fmin, const int fmax) const
{
    std::map<int,int>::const_iterator chIIter = fChannels.begin();
    std::map<int,int>::const_iterator chKIter = fChannels.begin();
    while(chKIter!=fChannels.end()) {
        o<<"\t"<<chKIter->first<<"                          ";
        chKIter++;
    }
    o<<std::endl;
 
    while(chIIter != fChannels.end()) {
        o<<chIIter->first<<"\t";
        std::map<int,int>::const_iterator chJIter = fChannels.begin(); 
        while(chJIter != fChannels.end()) {
            double corr=GetMeanCorrelation(chIIter->first,chJIter->first,fmin,fmax);
            o<<corr<<" ";
            chJIter++;
        }
        o<<std::endl;
        chIIter++;
    }  
}
 
 
void QChannelCovariance::Channels2Indices(const int chI, const int chJ, int& i, int& j) const
{
    i = j = 0;
    std::map<int,int>::const_iterator chIIter = fChannels.find(chI);
    std::map<int,int>::const_iterator chJIter = fChannels.find(chJ);
    if(chIIter == fChannels.end()) {
        QError err = QERR_SIZE_NOT_MATCH;
        std::stringstream msg;
        msg<<"Given channel "<<chI<<" not allowed";
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
        return;
    }
    if(chJIter == fChannels.end()) {
        QError err = QERR_SIZE_NOT_MATCH;
        std::stringstream msg;
        msg<<"Given channel "<<chJ<<" not allowed";
        err.SetDescription(__FILE__,__LINE__,msg.str());
        DianaThrow(err);
        return;
    }
    i = chIIter->second;
    j = chJIter->second;
}