QVector.cc

Go to the documentation of this file.

/*
* 
* Class QVector
*
*/
 
#include "QVector.hh"
#include "QVectorI.hh"
#include "QError.hh"
#include "QMatrix.hh"
#include <math.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_statistics_double.h>
#include <gsl/gsl_sort_vector.h>
#include <TH1D.h>
#include <TGraph.h>
#include <TPad.h>
#include <sstream>
#include <vector>
#include <exception>
#include <TMath.h>
 
#include "QRealComplexFFT.hh"
#ifdef _HAVE_FFTW3_
#include <fftw3.h>
#endif
 
QObjectImp(Diana::QVector);
 
Q_BEGIN_NAMESPACE
 
double*  QVector::ArrayAlloc(const UInt_t size) 
{
  double* a = 0;
  if(size > 0) {
    try{
#ifdef _HAVE_FFTW3A_
      a=(double *)fftw_malloc(size*sizeof(double));
#else
      a = new double[size];
#endif
    }
    catch(std::exception &e){
      std::stringstream msg;
      msg<<"QVector::ArrayAlloc: Attempt to allocate double array of size "
         << size << " failed!" << std::endl;
      msg<<"\n\tMessage: " << e.what();
      a=NULL;
    }
    if(a==NULL){
      std::stringstream msg;
      msg<<"QVector::ArrayAlloc: Attempt to allocate double array of size "
         << size << " failed!" << std::endl;
      DianaThrow(QError(QERR_UNDERFLOW,__FILE__,__LINE__,msg.str().c_str()));
    }
  }
  return a;
}
 
void QVector::ArrayCopy(const double* orig,const UInt_t origStride, double* dest,const UInt_t destStride,const UInt_t size)
{
  if(size > 0 && orig && dest) {
    if(origStride == 1 && destStride == 1)
      memcpy(dest,orig,sizeof(double)*size);
    else
      for(UInt_t i =0; i < size; i++) {
        dest[i*destStride] = orig[i*origStride]; 
      }
  }
}
 
void QVector::ArrayFree(double* array) 
{
  if(array){
#ifdef _HAVE_FFTW3A_
    fftw_free(array);
#else
    delete [] array;
#endif
  }
  array=NULL;
}
 
QVector::QVector()
  : fDataOwner(true),
    fSize(0),
    fAllocSize(0),
    fStride(1),
    fData(0),
    fMathVec(0)
{
  fMathVec=NULL;
}
 
QVector::QVector(const UInt_t size)
  : fDataOwner(true),
    fSize(0),
    fAllocSize(0),
    fStride(1),
    fData(0),
    fMathVec(0)
 
{
  Resize(size);
  fMathVec=NULL;
}
 
 
QVector::QVector(const QVector& in)
  : fDataOwner(true),
    fSize(0),
    fAllocSize(0),
    fStride(1),
    fData(0),
    fMathVec(0)
 
{
  Resize(in.Size());
  ArrayCopy(in.fData,in.fStride,fData,fStride,Size());
  fMathVec=NULL;
}
 
QVector::QVector(const QVector& orig,const bool dataOwner)
  : fDataOwner(dataOwner),
    fSize(0),
    fAllocSize(0),
    fStride(1),
    fData(0),
    fMathVec(0)
 
{
  if(dataOwner){
    Resize(orig.Size());
    ArrayCopy(orig.fData,orig.fStride,fData,fStride,Size());
  }
  else{
    fData=orig.fData;
    fStride=orig.fStride;
    fAllocSize=orig.fAllocSize;
    fSize=orig.fSize;
  }
  fMathVec=NULL;
}
 
QVector::QVector(double* vec,const UInt_t size,const UInt_t stride)
  : fDataOwner(false),
    fSize(size),
    fAllocSize(size),
    fStride(stride),
    fData(vec),
    fMathVec(0)
 
{
  fMathVec = NULL;
}
QVector::~QVector()
{
  if(fDataOwner && fData) {
    if(fData) ArrayFree(fData);
  }
  if(fMathVec) {
    free(fMathVec);
  }
  fData=NULL;
  fMathVec=NULL;
  fSize=0;
  fAllocSize=0;
}
 
void QVector::Resize(const UInt_t newSize){
  if(fDataOwner) {
    if(newSize > fAllocSize) {
      if(fData) ArrayFree(fData);
      fData = ArrayAlloc(newSize);
      fAllocSize = newSize;
      fStride=1;//check
    }
    fSize = newSize;
  } else if(newSize != (UInt_t)fSize) {
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                       "Not data owner, are you trying to resize a "
                       "QVector get from a QVector(Const)View?") );
  }
  if(fMathVec) free(fMathVec);
  fMathVec = NULL;
}
 
const QVector& QVector::operator=(const QVector& orig) 
{
    Resize(orig.fSize);
    ArrayCopy(orig.fData,orig.fStride,fData,fStride,fSize);    
    return *this;
}
const QVector& QVector::operator=(const QVectorI& orig) 
{
    Resize(orig.Size());
    for(UInt_t i = 0; i < fSize; i++) fData[i*fStride] = orig[i];
    return *this;
}
 
const QVector& QVector::operator=(const std::vector<double>& orig) 
{
    Resize(orig.size());
    for(UInt_t i = 0; i < fSize; i++) fData[i*fStride] = orig[i];
    return *this;
}
void QVector::SetMathVector() const
{
  if(!fMathVec) {
    fMathVec = (gsl_vector*) malloc(sizeof(gsl_vector));
    fMathVec->size = fSize;
    fMathVec->stride = fStride;
    fMathVec->data = fData;
    fMathVec->block = 0;
    fMathVec->owner = 0;
  }
}
void QVector::Initialize(const double val)
{
  for(UInt_t i = 0; i < fSize; i++) fData[i*fStride] = val;
}
 
double & QVector::operator[](const UInt_t i)
{
  if(i >= fSize) {
    std::stringstream stream;
    stream<<"Size: "<<fSize<<" Index: "<<i;
    DianaThrow( QError(QERR_OUT_OF_RANGE,__FILE__,__LINE__,stream.str()) );
  }
  return fData[i*fStride];
}
 
const double &QVector::operator[](const UInt_t i) const
{
  if(i >= fSize) {
    std::stringstream stream;
    stream<<"Size: "<<fSize<<" Index: "<<i;
    DianaThrow( QError(QERR_OUT_OF_RANGE,__FILE__,__LINE__,stream.str()) );
  }
  return fData[i*fStride];
}
 
 
QVector QVector::operator-() const 
{
    QVector tmp(*this);
    tmp*=(-1.);
    return tmp;
}
QVector QVector::Mult(const QVector&b) const
{
  QVector tmp(*this);
  b.SetMathVector();
  tmp.SetMathVector();
  gsl_vector_mul(tmp.fMathVec, b.fMathVec); 
  return tmp;
}
QVector QVector::Div(const QVector&b) const
{
  QVector tmp(*this);
  b.SetMathVector();
  tmp.SetMathVector();
  gsl_vector_div(tmp.fMathVec, b.fMathVec); 
  return tmp;
}
 
const QVector& QVector::operator/=(const double t) 
{
  return operator*=(1./t);
}
 
const QVector& QVector::operator*=(const double t) 
{
  SetMathVector();
  gsl_vector_scale(fMathVec,t);
  return *this;
}
 
 
const QVector & QVector::operator+=(const double v)
{
    SetMathVector();
    gsl_vector_add_constant(fMathVec,v);
    return *this;
}
 
const QVector & QVector::operator-=(const double v)
{
    SetMathVector();
    gsl_vector_add_constant(fMathVec,-v);
    return *this;
}
 
const QVector& QVector::operator+=(const QVector& other) 
{
    other.SetMathVector();
  SetMathVector();
  gsl_vector_add(fMathVec,other.fMathVec);
  return *this;
}
 
const QVector& QVector::operator-=(const QVector& other) 
{
  other.SetMathVector();
  SetMathVector();
  gsl_vector_sub(fMathVec,other.fMathVec);
  return *this;
}
 
const QVector & QVector::operator*=(const QVector &v) {
  v.SetMathVector();
  SetMathVector();
  gsl_vector_mul(fMathVec,v.fMathVec);
  return *this;
}
 
const QVector & QVector::operator/=(const QVector &v) {
  v.SetMathVector();
  SetMathVector();
  gsl_vector_div(fMathVec,v.fMathVec);
  return *this;
}
 
double QVector::operator*(const QVector& other) const
{
  other.SetMathVector();
  SetMathVector();
  double res;
  gsl_blas_ddot(fMathVec,other.fMathVec,&res);
  return res;
}
 
QVector QVector::operator+(const double t) const
{ 
  QVector tmp(*this);
  tmp+=t;
  return tmp;
}
 
QVector QVector::operator-(const double t) const
{
  QVector tmp(*this);
  tmp-=t;
  return tmp;
}
QVector QVector::operator*(const double t) const
{ 
  QVector tmp(*this);
  tmp*=t;
  return tmp;
}
 
QVector QVector::operator/(const double t) const
{
  QVector tmp(*this);
  tmp/=t;
  return tmp;
}
 
QVector QVector::operator+(const QVector& other) const
{
  QVector tmp(*this);
  tmp+=other;
  return tmp;
}
 
QVector QVector::operator-(const QVector& other) const
{
  QVector tmp(*this);
  tmp-=other;
  return tmp;
}
 
const QVector& QVector::Reverse()
{
        QVector output(Size());
        for(size_t i = 0; i < Size(); i++) {
                output[i] = operator[](Size()-1-i);
        }
        *this = output;
        return *this;
}
const QVector& QVector::Rebin(const int rebin,bool sum) {
  int newSize = Size()/rebin;
  double norm = (sum ? 1. : 1./rebin);
  for(int i=0;i<newSize;i++){
    (*this)[i]=norm*(*this)[rebin*i];
    for(int j=1;j<rebin;j++){
      (*this)[i]+=norm*(*this)[rebin*i+j];
    }
  }
  Resize(newSize);
  return *this;
}
const QVector& QVector::Shift(const int nstep)
{
        bool gtz = nstep > 0 ? true : false;
        if(nstep == 0 || Size() == 0) return *this;
        unsigned int unstep = abs(nstep) % Size();
        QVector rest(unstep);
  if(gtz) {
    for(unsigned int i = 0; i < unstep; i++) rest[i] = (*this)[Size() - unstep + i];
    for(int i = Size() - unstep - 1; i >= 0 ; i--) (*this)[i + unstep] = (*this)[i];
    for(unsigned int i = 0; i < unstep; i++) (*this)[i] = rest[i];
  } else {      
    for(unsigned int i = 0; i < unstep; i++) rest[i] = (*this)[i];
    for(unsigned int i = unstep; i < Size(); i++) (*this)[i - unstep] = (*this)[i];
    for(unsigned int i = 0; i < unstep; i++) (*this)[i + Size() - unstep] = rest[i];
  }
  return *this;
}
 
const QVector& QVector::ShiftReal(const double fstep)
{
  int nstep = (int) ceil(fstep);
  Shift(nstep);
  double fracs = nstep - fstep;
  double zero = (*this)[0];
  for(int i = 0; i < (int)Size() - 1; i++) {
    (*this)[i] = (*this)[i]*(1.-fracs) + fracs*(*this)[i+1]; 
  }
  (*this)[Size()-1] = (*this)[Size()-1]*(1.-fracs) + fracs*zero;
  
  return *this;
}
 
QVector QVector::Abs() const {
  QVector vec(*this);
  for(UInt_t i=0;i<vec.Size();i++)
    vec[i]=TMath::Abs(vec[i]);
  return vec;
}
 
const QVector& QVector::Sort(const bool asc)
{
  SetMathVector();
  gsl_sort_vector(fMathVec);
  if(!asc) gsl_vector_reverse(fMathVec);
  return *this;
}
 
const QVector& QVector::Differentiate()
{
  for(UInt_t i = Size()-1; i > 0; i--)
    (*this)[i] = (*this)[i] - (*this)[i-1] ;    
  (*this)[0] = (*this)[1];
  
  return *this;
}
 
QVector QVector::Derivative2P(const double delta) const
{
  QVector deriv = (*this)/delta;
  deriv.Differentiate();
  deriv[0] = deriv[1]; 
  return deriv;
}
 
QVector QVector::Derivative2PF(const double delta) const
{
  const QVector &me = (*this);
  QVector result(me.Size());
  for(UInt_t i = 0; i < me.Size()-1; i ++) {
    result[i] = (me[i+1]-me[i])/delta;
  }
  result[me.Size()-1] = result[me.Size()-2];
  return result;
}
QVector QVector::Derivative3P(const double delta) const
{
  const QVector &me = (*this);
  QVector result(me.Size());
  result[0] = (-3.*me[0]+4.*me[1]-me[2]) / (2.*delta);
  for(UInt_t i = 1; i < me.Size()-1; i ++) {
    result[i] = (me[i+1]-me[i-1]) / (2.*delta);
  }
  result[me.Size()-1] = (me[me.Size()-3]-4.*me[me.Size()-2]+3.*me[me.Size()-1]) / (2.*delta);
  return result;
}
 
QVector QVector::Derivative3PB(const double delta) const
{
  const QVector &me = (*this);
  QVector result(me.Size());
  result[0] = (-3.*me[0]+4.*me[1]-me[2]) / (2.*delta);
  result[1] = (me[2]-me[0]) / (2.*delta);
  for(UInt_t i = 2; i < me.Size(); i ++) {
    result[i] = (me[i-2]-4.*me[i-1]+3.*me[i]) / (2.*delta);
  }
  return result;
}
 
QVector QVector::Derivative3PF(const double delta) const
{
  const QVector &me = (*this);
  QVector result(me.Size());
  UInt_t i = 0;
  for(i = 0; i < me.Size()-2; i++) {
    result[i] = (-3.*me[i]+4.*me[i+1]-me[i+2]) / (2.*delta);
  }
  i = me.Size()-2;
  result[i] = (me[i+1]-me[i-1]) / (2.*delta);
  i = me.Size()-1;
  result[i] = (me[i-2]-4.*me[i-1]+3.*me[i]) / (2.*delta);
  return result;
}
 
QVector QVector::Derivative5P(const double delta) const
{
  const QVector &me = (*this);
  const UInt_t size = me.Size();
  QVector result(size);
  result[0] = (-25.*me[0]+48.*me[1]-36.*me[2]+16.*me[3]-3.*me[4]) / (12.*delta);
  result[1] = (-3.*me[0]-10.*me[1]+18.*me[2]-6.*me[3]+1.*me[4]) / (12.*delta);
  for(UInt_t i = 2; i < size-2; i++) {
    result[i] = (me[i-2]-8.*me[i-1]+8.*me[i+1]-me[i+2]) / (12.*delta);
  }
  result[size-2] = (3.*me[size-1]+10.*me[size-2]-18.*me[size-3]+6.*me[size-4]-1.*me[size-5]) / (12.*delta);
  result[size-1] = (25.*me[size-1]-48.*me[size-2]+36.*me[size-3]-16.*me[size-4]+3.*me[size-5]) / (12.*delta);
  
  return result;
}
 
QVector QVector::Derivative5PB(const double delta) const
{
  const QVector &me = (*this);
  const UInt_t size = me.Size();
  QVector result(size);
  result[0] = (-25.*me[0]+48.*me[1]-36.*me[2]+16.*me[3]-3.*me[4]) / (12.*delta);
  result[1] = (-3.*me[0]-10.*me[1]+18.*me[2]-6.*me[3]+1.*me[4]) / (12.*delta);
  for(UInt_t i = 2; i < 4; i++) {
    result[i] = (me[i-2]-8.*me[i-1]+8.*me[i+1]-me[i+2]) / (12.*delta);
    }
  for(UInt_t i = 4; i < size; i++)
    result[i] = (25.*me[i]-48.*me[i-1]+36.*me[i-2]-16.*me[i-3]+3.*me[i-4]) / (12.*delta);
  
  return result;
}
 
QVector QVector::DerivativeEulero(const double delta) const
{
  
  const QVector &me = (*this);
  const UInt_t size = me.Size();
  QVector result(size);
    result[0] = (-3.*me[0]+4.*me[1]-me[2]) / (2.*delta);
    for(UInt_t i = 1; i < size; i++) 
      result[i] = (me[i]-me[i-1])*2./delta - result[i-1];
    return result;
}
 
 
const QVector& QVector::Integrate(const double constTerm)
{
  (*this)[0] += constTerm;
  for(UInt_t i = 1; i < Size(); i++)
    (*this)[i] += (*this)[i-1];
  return *this;
}
 
QMatrix QVector::T() const
{
  QMatrix tmp(1,Size());
  tmp.SetRow(0,*this);
  return tmp;
}
 
QMatrix QVector::Reshape(UInt_t n,UInt_t m) {
 
  //  IMPORANT NOTE: These algorithms do not work if the vector has a non-unity stride!
 
  if(n*m != Size())
    throw QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                 Form("Cannot reshape vector of size %ui into matrix of size %uix%ui!",Size(),n,m));
  if(fStride != 1)
    throw QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                 "Cannot reshape vector with a stride greater than 1!");
  QMatrix Ret(fData,n,m);
  return Ret;
}
 
const QMatrix QVector::Reshape(UInt_t n,UInt_t m) const {
 
//  IMPORANT NOTE: These algorithms do not work if the vector has a non-unity stride!
 
  if(n*m != Size())
    throw QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                 Form("Cannot reshape vector of size %ui into matrix of size %uix%ui!",Size(),n,m));
  if(fStride != 1)
    throw QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                 "Cannot reshape vector with a stride greater than 1!");
  QMatrix Ret(fData,n,m);
  return Ret;
}
 
 
void QVector::Append(const double val)
{
  if(fDataOwner) {
    double* temp = ArrayAlloc(Size()+1);
    ArrayCopy(fData,fStride,temp,1,Size());
    if(fData) ArrayFree(fData);
    //  fData = ArrayAlloc(newSize);
    fSize=Size()+1;
    fAllocSize = fSize;;
    fStride=1;//check
    temp[Size()-1]=val;
    fData=temp;
    fStride=1;
    if(fMathVec) free(fMathVec);
    fMathVec = NULL;
  } 
  else {
    DianaThrow( QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,"Not data owner, you are trying to resize a QVector using Append") );
  }
}
 
double QVector::Norma() const
{
  return sqrt((*this)*(*this));
}
 
 
double QVector::Sum(const UInt_t nelem,const UInt_t first) const
{
  double sum = 0;
  UInt_t last = first + nelem;
  for(UInt_t i = first; i < last; i++) sum += this->operator[](i);
  return sum;
}
 
double QVector::GetRMS(const UInt_t nelem,const UInt_t first) const
{
  SetMathVector();
  double mean = Sum(nelem, first);
  mean /= nelem;
  double rms = gsl_stats_variance_m(&(fMathVec->data[first]), 1,nelem, mean);
  return sqrt(rms);
}
 
double QVector::GetMedian() const
{
  
  SetMathVector();
  double median = 0;
  if (Size() > 0) {
    // make a copy of the gsl vector and sort the copy in place
    gsl_vector* copy = gsl_vector_alloc(Size());
    gsl_vector_memcpy(copy, fMathVec);
    gsl_sort_vector(copy);
    
    if (Size() % 2 == 1) { // odd number of elements
      median = gsl_vector_get(copy, (Size() - 1) / 2);
    }
    else { // even number of elements
      median = (  gsl_vector_get(copy, Size() / 2 - 1)
                  + gsl_vector_get(copy, Size() / 2    ) ) / 2.0;
    }
    
    // free the memory for the copy
    gsl_vector_free(copy);
  }
  return median;
}
 
double QVector::GetMedianAbsoluteDeviation() const
{
  double mad = -1;
  if (Size() > 0) {
    double median = GetMedian();
    QVector absDeviations(Size());
    for (unsigned int i = 0; i < Size(); ++i) {
      absDeviations[i] = fabs(operator[](i) - median);
    }
    mad = absDeviations.GetMedian();
  }
  return mad;
}
 
int QVector::GetMaxIndex(const UInt_t nelem,const UInt_t first) const
{
  
  SetMathVector();
  gsl_vector_const_view vectorview = gsl_vector_const_subvector(fMathVec,first,nelem);
  return gsl_vector_max_index(&(vectorview.vector))+first;
}
 
int QVector::GetMinIndex(const UInt_t nelem,const UInt_t first) const
{
  
  SetMathVector();
  gsl_vector_const_view vectorview = gsl_vector_const_subvector(fMathVec,first,nelem);
  return gsl_vector_min_index(&(vectorview.vector))+first;
}
 
double* QVector::GetArray() const {
  double *a = ArrayAlloc(fSize);
  ArrayCopy(fData,fStride,a,1,fSize);
  return a;
}
QVector QVector::GetSubVector(UInt_t size,UInt_t start,UInt_t stride) {
  if(start+size > Size()) 
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Sub vector would exceed beyond original vector!"));
  stride = (stride < 1 ? 1 : stride);
  return QVector(fData+fStride*start,size,fStride*stride);
}
 
const QVector QVector::GetSubVector(UInt_t size,UInt_t start,UInt_t stride)
  const {
  if(start+size > Size()) 
    DianaThrow(QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                      "Sub vector would exceed beyond original vector!"));
  stride = (stride < 1 ? 1 : stride);
  return QVector(fData+fStride*start,size,fStride*stride);
}
 
 
void QVector::SetArray(const double* orig,const UInt_t size) 
{
  Resize(size);
  ArrayCopy(orig,1,fData,fStride,size);
}
 
 
std::string QVector::sqlString() const{
  std::string s = "";
  char d[20];
  if(Size() > 0)  {
    s += "'{";
    for (UInt_t i=0;i<Size();i++){
      if(i!=0)s+=std::string(",");
      // retain 10 significant digis
      snprintf(d,20,"%.10g",operator[](i));
      s+=d;    
    }
    s+=std::string("}'");
  }
  else
    s+="NULL"; 
  return s;
}
 
void QVector::Draw(Option_t* option)
{
  TH1D* h =  new TH1D("Diana::QVector","",fSize,-0.5,double(fSize)-0.5);
  h->SetDirectory(0);
  for(UInt_t i = 0; i < fSize; i++)
    h->SetBinContent(i+1,this->operator[](i));
  h->SetBit(kCanDelete);
  h->Draw(option);
  
  if(gPad) {
    gPad->Update();
    gPad->Modified();
  }
  
}
 
 
TGraph* QVector::GetGraph(double samplingFrequency, double scale) const
{
  const QVector& timeline = *this;
  TGraph* graph = new TGraph(timeline.Size());
  for(UInt_t i = 0; i < timeline.Size(); i++) 
    {
      graph->SetPoint(i,double(i)/samplingFrequency,timeline[i]*scale);
    }
  return graph;
}
 
void QVector::Dump(std::ostream&s) const
{
  s<<"Size: "<<this->Size()<<"; Elements: ";
  for (UInt_t i = 0; i < this->Size(); i++) s<<this->operator[](i)<<" ";
  s<<std::endl;
}
 
 
Q_END_NAMESPACE
 
 
/********************************************
  NON MEMBER FUNCTIONS
*********************************************/
 
Diana::QVector operator+(double t, const Diana::QVector &v){
  Diana::QVector tmp(v);
  tmp+=t;
  return tmp;
}
 
Diana::QVector operator-(double t, const Diana::QVector &v){
  Diana::QVector tmp(v);
  tmp-=t;
  return -tmp;
}
Diana::QVector operator*(double t, const Diana::QVector &v){
  Diana::QVector tmp(v);
  tmp*=t;
  return tmp;
}
 
Diana::QVector operator/(double t, const Diana::QVector &v){
  
  Diana::QVector tmp(v.Size());
  tmp.Initialize(t);
  return tmp.Div(v);
}
 
Diana::QVector abs(const Diana::QVector &v){
  return v.Abs();
}
 
Diana::QVector sin(const Diana::QVector &v) {
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i]=sin(v[i]);
  return tmp;
}
 
Diana::QVector cos(const Diana::QVector &v) {
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i]=cos(v[i]);
  return tmp;
}
 
Diana::QVector tan(const Diana::QVector &v) {
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i]=tan(v[i]);
  return tmp;
}
 
Diana::QVector exp(const Diana::QVector &v) {
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i]=exp(v[i]);
  return tmp;
}
 
Diana::QVector log(const Diana::QVector &v) {
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i] = TMath::Log(v[i]);
  return tmp;
}
 
Diana::QVector log10(const Diana::QVector &v) {
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i] = TMath::Log10(v[i]);
  return tmp;
}
 
 
Diana::QVector pow(const Diana::QVector &v,double p) {
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i]=TMath::Power(v[i],p);
  return tmp;
}
 
Diana::QVector pow(double p,const Diana::QVector &v) {
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i]=TMath::Power(p,v[i]);
  return tmp;
}
 
Diana::QVector pow(const Diana::QVector &u,const Diana::QVector &v) {
  if(u.Size()!=v.Size())
    throw QError(QERR_SIZE_NOT_MATCH,__FILE__,__LINE__,
                 "Cannot raise a vector to the power of another vector if the sizes do not match!");
  Diana::QVector tmp(v.Size());
  for(UInt_t i=0;i<v.Size();i++) tmp[i]=TMath::Power(u[i],v[i]);
  return tmp;
}
 
 
double max(const Diana::QVector &v){
  return v.GetMax();
}
 
double min(const Diana::QVector &v){
  return v.GetMin();
}