QVectorTools.cc

Go to the documentation of this file.

#include "QVectorTools.hh"
 
#include "QRealComplexFFT.hh"
#include <TRandom.h>
#include <TMath.h>
#include <math.h>
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::Zeros(UInt_t size) {
  Diana::QVector vec(size);
  vec.Initialize(0);
  return vec;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::Range(UInt_t size) {
  Diana::QVector vec(size);
  for(UInt_t i=0;i<size;i++) vec[i]=i;
  return vec;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::Range(UInt_t min,UInt_t max,UInt_t stride) {
  UInt_t NP=(max-min-1)/stride+1;
  Diana::QVector vec(NP);
  for(UInt_t i=0;i<NP;i++) vec[i]=min+i*stride;
  return vec;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::LinSpace(Double_t min,Double_t max,UInt_t size) {
  return (max-min)*Range(size)/(size-1)+min;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::LogSpace(Double_t min,Double_t max,UInt_t size) {
  return pow(10.,LinSpace(min,max,size));
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::Random::Gauss(UInt_t Size,double mean,double std) {
  Diana::QVector t(Size);
  for(UInt_t i=0;i<Size;i++)
    t[i] = gRandom->Gaus(mean,std);
  return t;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::Random::Uniform(UInt_t Size,double max) {
  return Uniform(Size,0,max);
}
    
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::Random::Uniform(UInt_t Size,double min,double max) {
  Diana::QVector t(Size);
  for(UInt_t i=0;i<Size;i++)
    t[i] = gRandom->Uniform(min,max);
  return t;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::Random::NoiseWaveformFromANPS(const Diana::QVector &anps) {
 
  const double sqrt2over2 = sqrt(0.5);
  
  Diana::QVector noise;
  Diana::QVectorC noiseC(anps.Size());
 
  noiseC[0].Set(0,0); // Set zero DC value
  noiseC[noiseC.Size()-1] = sqrt2over2*gRandom->Exp(sqrt(anps[noiseC.Size()-1]));
  for(size_t i=1;i<noiseC.Size()-1;i++)
    noiseC[i].SetPolar(sqrt2over2*gRandom->Exp(sqrt(anps[i])),
                       TMath::TwoPi()*gRandom->Uniform(0,1));
 
  Diana::gRealComplexFFT->TransformFromFreq(noiseC,noise);
 
  return noise;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::Convolve(const Diana::QVector &vecA,const Diana::QVector &vecB,ConvolveReturn mode) {
  Diana::QVector output = Zeros(vecA.Size()+vecB.Size()-1);
  for(size_t i=0;i<output.Size();i++) {
    for(size_t j=0;j<vecA.Size();j++) {
      if(((int)i-(int)j)<0 || (i-j)>=vecB.Size()) continue;
      output[i] += vecA[j]*vecB[i-j];
    }
  }
  if(mode == kConvolveSame) {
    int retSize = std::max(vecA.Size(),vecB.Size());
    int offset = (output.Size()-retSize)/2;
    output=output.GetSubVector(retSize,offset);
  }
  else if(mode == kConvolveValid) {
    int retSize = std::max(vecA.Size(),vecB.Size())-std::min(vecA.Size(),vecB.Size())+1;
    int offset = (output.Size()-retSize)/2;
    output=output.GetSubVector(retSize,offset);
  }
  return output;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::CircularConvolve(const Diana::QVector &vecA,const Diana::QVector &vecB) {
  Diana::QVector output = Zeros(vecA.Size()+vecB.Size()-1);
  const Diana::QVector &longVec = (vecA.Size() >= vecB.Size() ? vecA : vecB);
  const Diana::QVector &shortVec = (vecA.Size() < vecB.Size() ? vecA : vecB);
  for(size_t i=0;i<output.Size();i++) {
    for(size_t j=0;j<shortVec.Size();j++) {
      output[i] += shortVec[j]*longVec[(i-j+longVec.Size())%longVec.Size()];
    }
  }
  return output;
}
 
//_____________________________________________________________________________
Diana::QVector Diana::QVectorTools::ConvolveFFT(const Diana::QVector &vecA,const Diana::QVector &vecB,ConvolveReturn mode) {
  
  const Int_t SizeToOverlapAdd = 5000;
  
  Diana::QVectorC vecA_f,vecB_f;
  Diana::QRealComplexFFT FFT;
  Diana::QVector output;
  if(std::abs((int)(vecA.Size()-vecB.Size()))<SizeToOverlapAdd) {
    // If the sizes are very close, just pad them, fft them and add them.
    int padSize = vecA.Size()+vecB.Size()-1;
    Diana::QVector vecA_p = Diana::QFFT::ZeroPad(vecA,padSize-vecA.Size(),Diana::QFFT::kRight);
    Diana::QVector vecB_p = Diana::QFFT::ZeroPad(vecB,padSize-vecB.Size(),Diana::QFFT::kRight);
    FFT.TransformToFreq(vecA_p,vecA_f);
    FFT.TransformToFreq(vecB_p,vecB_f);
    Diana::QVectorC output_f = vecA_f.Mult(vecB_f);
    FFT.TransformFromFreq(output_f,output);
    output=output.GetSubVector(vecA.Size()+vecB.Size()-1);
  }
  else {
    // if the sizes are vastly different, use the overlap and add method
    output.Resize(vecA.Size()+vecB.Size()-1);
    output.Initialize(0);
 
    const Diana::QVector &longVec  = (vecA.Size() > vecB.Size() ? vecA : vecB);
    const Diana::QVector &shortVec = (vecA.Size() < vecB.Size() ? vecA : vecB);
 
    UInt_t longSize = longVec.Size();
    UInt_t shortSize = shortVec.Size();
    UInt_t skip,dum;
    Diana::QFFT::IsPowerOf2(shortSize,skip,dum);
    skip=pow(2,skip);
 
    Diana::QVector output_seg;
    Diana::QVector shortVec_p = Diana::QFFT::ZeroPad(shortVec,skip-1,Diana::QFFT::kRight);
    Diana::QVectorC shortVec_f,long_seg_f,output_seg_f;
    FFT.TransformToFreq(shortVec_p,shortVec_f);
 
    for(UInt_t segStart=0;segStart<longSize;segStart+=skip) {
      if(segStart+skip<longSize) {
        Diana::QVector long_seg_p = Diana::QFFT::ZeroPad(longVec.GetSubVector(skip,segStart),shortSize-1,Diana::QFFT::kRight);
        FFT.TransformToFreq(long_seg_p,long_seg_f);
        output_seg_f = long_seg_f.Mult(shortVec_f);
        FFT.TransformFromFreq(output_seg_f,output_seg);
        output.GetSubVector(shortSize+skip-1,segStart) += output_seg;
      }
      else {
        output_seg=ConvolveFFT(shortVec,longVec.GetSubVector(longSize-segStart,segStart));
        output.GetSubVector(output.Size()-segStart,segStart) += output_seg;
      }
    }
  }
  if(mode == kConvolveSame) {
    int retSize = std::max(vecA.Size(),vecB.Size());
    int offset = (output.Size()-retSize)/2;
    output=output.GetSubVector(retSize,offset);
  }
  else if(mode == kConvolveValid) {
    int retSize = std::max(vecA.Size(),vecB.Size())-std::min(vecA.Size(),vecB.Size())+1;
    int offset = (output.Size()-retSize)/2;
    output=output.GetSubVector(retSize,offset);
  }  
 
  return output;
}