NcTransform.cxx

Go to the documentation of this file.

 
#include "NcTransform.h"
#include "Riostream.h"
 
ClassImp(NcTransform); // Class implementation to enable ROOT I/O
 
NcTransform::NcTransform(const char* name,const char* title) : TNamed(name,title)
{
 
 Reset();
 fSample=0;
}

NcTransform::~NcTransform()
{
 
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
}

NcTransform::NcTransform(const NcTransform& q) : TNamed(q)
{
 
 fProc=0;
 fN=q.fN;
 fReIn=q.fReIn;
 fImIn=q.fImIn;
 fReOut=q.fReOut;
 fImOut=q.fImOut;
 fSample=q.fSample;
}

void NcTransform::Reset()
{
 
 fProc=0;
 fN=0;
 fReIn.Set(0);
 fImIn.Set(0);
 fReOut.Set(0);
 fImOut.Set(0);
}

void NcTransform::SetSamplingFrequency(Float_t f)
{
 
 fSample=f;
}

Float_t NcTransform::GetSamplingFrequency() const
{
 
 return fSample;
}

void NcTransform::Load(Int_t n,Double_t* re,Double_t* im,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (n<1) return;
 
 fN=n;
 if (re) fReIn.Set(n);
 if (im) fImIn.Set(n);
 
 for (Int_t i=0; i<n; i++)
 {
  if (re) fReIn.SetAt(re[i],i);
  if (im) fImIn.SetAt(im[i],i);
 }
}

void NcTransform::Load(TArray* re,TArray* im,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 Int_t nre=0;
 Int_t nim=0;
 if (re) nre=re->GetSize();
 if (im) nim=im->GetSize();
 
 Int_t n=nre;
 if (!n) n=nim;
 if (nre && nim>nre) n=nre;
 if (nim && nre>nim) n=nim;
 
 if (n<1) return;
 
 fN=n;
 if (nre) fReIn.Set(n);
 if (nim) fImIn.Set(n);
 
 Double_t valre=0;
 Double_t valim=0;
 for (Int_t i=0; i<n; i++)
 {
  if (nre)
  {
   valre=re->GetAt(i);
   fReIn.SetAt(valre,i);
  }
  if (nim)
  {
   valim=im->GetAt(i);
   fImIn.SetAt(valim,i);
  }
 }
}

void NcTransform::Load(NcSample* s,Int_t i,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (!s)
 {
  cout << " *" << ClassName() <<"::Load* No input sample specified." << endl;
  return;
 }
 
 Int_t n=s->GetN();
 Int_t store=s->GetStoreMode();
 Int_t dim=s->GetDimension();
 
 if (n<1 || !store || dim<1 || i<1 || i>dim)
 {
  cout << " *" << ClassName() <<"::Load* Inconsistent input for NcSample treatment." << endl;
  cout << " Store Mode:" << store << " Entries:" << n << " Dimension:" << dim << " i:" << i << " f:" << fSample << endl; 
  return;
 }
 
 fN=n;
 fReIn.Set(n);
 
 Double_t val=0;
 for (Int_t idx=1; idx<=n; idx++)
 {
  val=s->GetEntry(idx,i);
  fReIn.SetAt(val,idx-1);
 }
}

void NcTransform::Load(NcSample* s,TString name,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (!s)
 {
  cout << " *" << ClassName() <<"::Load* No input sample specified." << endl;
  return;
 }
 
 Int_t n=s->GetN();
 Int_t store=s->GetStoreMode();
 Int_t dim=s->GetDimension();
 Int_t i=s->GetIndex(name);
 
 
 if (n<1 || !store || dim<1 || !i)
 {
  cout << " *" << ClassName() <<"::Load* Inconsistent input for NcSample treatment." << endl;
  cout << " Store Mode:" << store << " Entries:" << n << " Dimension:" << dim << " name:" << name << " f:" << fSample << endl; 
  return;
 }
 
 Load(s,i,f);
}

void NcTransform::Load(TH1* h,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (!h)
 {
  cout << " *" << ClassName() <<"::Load* No input histogram specified." << endl;
  return;
 }
 
 Int_t nbins=h->GetNbinsX();
 Double_t nentries=h->GetEntries();
 
 if (!nbins || !nentries)
 {
  cout << " *" << ClassName() <<"::Load* Inconsistent input for histogram treatment." << endl;
  cout << " Nbins:" << nbins << " Nentries:" << nentries << " f:" << fSample << endl; 
  return;
 }
 
 fN=nbins;
 fReIn.Set(nbins);
 
 Double_t val=0;
 for (Int_t i=1; i<=nbins; i++)
 {
  val=h->GetBinContent(i);
  fReIn.SetAt(val,i-1);
 }
}

void NcTransform::Load(TGraph* gr,Float_t f)
{
 
 Reset();
 
 if (f>=0) fSample=f;
 
 if (!gr)
 {
  cout << " *" << ClassName() <<"::Load* No input TGraph object specified." << endl;
  return;
 }
 
 Int_t n=gr->GetN();
 
 if (!n)
 {
  cout << " *" << ClassName() <<"::Load* Inconsistent input for TGraph treatment." << endl;
  cout << " n:" << n << " f:" << fSample << endl; 
  return;
 }
 
 fN=n;
 fReIn.Set(n);
 
 gr->Sort();
 
 Double_t x=0;
 Double_t y=0;
 for (Int_t i=0; i<n; i++)
 {
  gr->GetPoint(i,x,y);
  fReIn.SetAt(y,i);
 }
}

void NcTransform::LoadResult()
{
 
 fReIn=fReOut;
 fImIn=fImOut;
 fReOut.Set(0);
 fImOut.Set(0);
}

Int_t NcTransform::GetN() const
{
 
 return fN;
}

TArrayD NcTransform::GetData(TString sel) const
{
 
 if (sel.Contains("RE") && sel.Contains("in")) return fReIn;
 if (sel.Contains("IM") && sel.Contains("in")) return fImIn;
 if (sel.Contains("RE") && sel.Contains("out")) return fReOut;
 if (sel.Contains("IM") && sel.Contains("out")) return fImOut;
 
 TArrayD data(fN);
 Double_t pi=acos(-1.);
 Double_t re=0;
 Double_t im=0;
 Double_t amp=0;
 Double_t phi=0;
 for (Int_t i=0; i<fN; i++)
 {
  if (sel.Contains("in"))
  {
   re=fReIn.At(i);
   im=fImIn.At(i);
  }
  if (sel.Contains("out"))
  {
   re=fReOut.At(i);
   im=fImOut.At(i);
  }
  amp=sqrt(re*re+im*im);
  phi=0;
  if (im || re) phi=atan2(im,re);
 
  if (sel.Contains("AMP")) data.SetAt(amp,i);
  if (sel.Contains("PHIR")) data.SetAt(phi,i);
  if (sel.Contains("PHID")) data.SetAt(phi*180./pi,i);
 }
 
 return data;
}

void NcTransform::Fourier(TString mode,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 
 if (fN<1) return;
 
 Int_t n=1+fN/2;
 Float_t maxfrac=0.5;
 if (sel.Contains("n") || sel.Contains("t") || sel.Contains("2"))
 {
  n=fN;
  maxfrac=1;
 }
 
 // Construct the Fast Fourier Transform (FFT) processor
 TString opt=mode;
 // Comply with the TVirtualFFT conventions
 if (mode=="C2C") opt="C2CFORWARD";
 if (mode=="C2CI") opt="C2CBACKWARD";
 // Set the optimalisation initialisation to "estimate" and create a new processor
 opt+=" ES K";
 
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 
 fProc=TVirtualFFT::FFT(1,&fN,opt);
 
 // Enter the input data
 Int_t nReIn=fReIn.GetSize();
 Int_t nImIn=fImIn.GetSize();
 if (mode=="R2C")
 {
  Double_t* data=fReIn.GetArray();
  fProc->SetPoints(data);
 }
 else
 {
  for (Int_t i=0; i<fN; i++)
  {
   if (nReIn && nImIn) fProc->SetPoint(i,fReIn[i],fImIn[i]);
   if (nReIn && !nImIn) fProc->SetPoint(i,fReIn[i],0);
   if (!nReIn && nImIn) fProc->SetPoint(i,0,fImIn[i]);
  }
 }
 
 // Perform the Fast Fourier Transform
 fProc->Transform();
 
 Double_t rN=fN;
 
 // Copy the resulting data in the output arrays
 fReOut.Set(fN);
 fImOut.Set(fN);
 Double_t re=0;
 Double_t im=0;
 for (Int_t i=0; i<fN; i++)
 {
  fProc->GetPointComplex(i,re,im);
  re/=sqrt(rN);
  im/=sqrt(rN);
  fReOut.SetAt(re,i);
  fImOut.SetAt(im,i);
 }
 
 if (!hist) return;
 
 if ((sel.Contains("Hz") || sel.Contains("t")) && fSample<=0) return;
 
 // Initialize the histogram title
 TString title="";
 if (mode=="C2R" || mode=="C2CI") title+="Inverse ";
 title+="DFT (";
 title+=mode;
 title+=") ";
 
 TString title2="";
 
 // Define and fill the requested result histogram
 if (sel.Contains("k"))
 {
  hist->SetBins(n,0,n-1);
  title+="index frequency domain";
  if (mode=="C2R" || mode=="C2CI")
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Index k";
  if (sel.Contains("RE")) title+=";Re(Q[k])";
  if (sel.Contains("IM")) title+=";Im(Q[k])";
  if (sel.Contains("AMP")) title+=";Amplitude |Q[k]|";
  if (sel.Contains("PHIR")) title+=";Phase #varphi[k] (rad)";
  if (sel.Contains("PHID")) title+=";Phase #varphi[k] (deg)";
 }
 if (sel.Contains("f"))
 {
  hist->SetBins(n,0,maxfrac);
  title+="fractional frequency domain";
  if (mode=="C2R" || mode=="C2CI")
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Fraction f of sampling rate";
  if (sel.Contains("RE")) title+=";Re(Q[f])";
  if (sel.Contains("IM")) title+=";Im(Q[f])";
  if (sel.Contains("AMP")) title+=";Amplitude |Q[f]|";
  if (sel.Contains("PHIR")) title+=";Phase #varphi[f] (rad)";
  if (sel.Contains("PHID")) title+=";Phase #varphi[f] (deg)";
 }
 if (sel.Contains("Hz"))
 {
  hist->SetBins(n,0,maxfrac*fSample);
  title+="actual frequency domain";
  if (mode=="C2R" || mode=="C2CI")
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Frequency #nu (Hz)";
  if (sel.Contains("RE")) title+=";Re(Q[#nu])";
  if (sel.Contains("IM")) title+=";Im(Q[#nu])";
  if (sel.Contains("AMP")) title+=";Amplitude |Q[#nu]|";
  if (sel.Contains("PHIR")) title+=";Phase #varphi[#nu] (rad)";
  if (sel.Contains("PHID")) title+=";Phase #varphi[#nu] (deg)";
 }
 
 if (sel.Contains("n"))
 {
  hist->SetBins(fN,0,fN);
  title+="sampling time domain";
  if (mode=="R2C" || mode=="C2C") title+=" input";
  if (fSample>0) title2.Form(" (%-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Sample number n";
  if (mode=="R2C" || mode=="C2R")
  {
   title+=";Value X[n]";
  }
  else
  {
   if (sel.Contains("RE")) title+=";Re(X[n])";
   if (sel.Contains("IM")) title+=";Im(X[n])";
   if (sel.Contains("AMP")) title+=";Amplitude |X[n]|";
   if (sel.Contains("PHIR")) title+=";Phase #varphi[n] (rad)";
   if (sel.Contains("PHID")) title+=";Phase #varphi[n] (deg)";
  }
 }
 if (sel.Contains("t"))
 {
  hist->SetBins(fN,0,rN/fSample);
  title+="actual time domain";
  if (mode=="R2C" || mode=="C2C") title+=" input";
  if (fSample>0) title2.Form(" (%-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Time t (seconds)";
  if (mode=="R2C" || mode=="C2R")
  {
   title+=";Value X[t]";
  }
  else
  {
   if (sel.Contains("RE")) title+=";Re(X[t])";
   if (sel.Contains("IM")) title+=";Im(X[t])";
   if (sel.Contains("AMP")) title+=";Amplitude |X[t]|";
   if (sel.Contains("PHIR")) title+=";Phase #varphi[t] (rad)";
   if (sel.Contains("PHID")) title+=";Phase #varphi[t] (deg)";
  }
 }
 
 hist->SetTitle(title);
 
 Double_t pi=acos(-1.);
 Double_t amp=0;
 Double_t phi=0;
 re=0;
 im=0;
 for (Int_t i=0; i<n; i++)
 {
  if (sel.Contains("n") || sel.Contains("t")) // Time domain data requested
  {
   if (mode=="R2C")
   {
    if (nReIn) hist->SetBinContent(i+1,fReIn.At(i));
    continue;
   }
   if (mode=="C2R")
   {
    hist->SetBinContent(i+1,fReOut.At(i));
    continue;
   }
   if (mode=="C2C") 
   {
    if (nReIn) re=fReIn.At(i);
    if (nImIn) im=fImIn.At(i);
   }
   if (mode=="C2CI") 
   {
    re=fReOut.At(i);
    im=fImOut.At(i);
   }
  }
  else // Frequency domain data requested
  {
   if (mode=="C2R" || mode=="C2CI") // Inverse transformation
   {
    if (nReIn) re=fReIn.At(i);
    if (nImIn) im=fImIn.At(i);
   }
   else // Forward transformation
   {
    re=fReOut.At(i);
    im=fImOut.At(i);
   }
  }
  amp=sqrt(re*re+im*im);
  phi=0;
  if (im || re) phi=atan2(im,re);
 
  if (sel.Contains("RE")) hist->SetBinContent(i+1,re);
  if (sel.Contains("IM")) hist->SetBinContent(i+1,im);
  if (sel.Contains("AMP")) hist->SetBinContent(i+1,amp);
  if (sel.Contains("PHIR")) hist->SetBinContent(i+1,phi);
  if (sel.Contains("PHID")) hist->SetBinContent(i+1,phi*180./pi);
 }
}

void NcTransform::Hartley(Int_t mode,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 
 if (!mode || fN<1) return;
 
 Int_t n=1+fN/2;
 Float_t maxfrac=0.5;
 if (sel.Contains("n") || sel.Contains("t") || sel.Contains("2"))
 {
  n=fN;
  maxfrac=1;
 }
 
 // Construct the Fast Fourier Transform (FFT) processor
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 
 fProc=TVirtualFFT::FFT(1,&fN,"DHT ES K");
 
 // Enter the input data
 Double_t* data=fReIn.GetArray();
 fProc->SetPoints(data);
 
 // Perform the Discrete Hartley Transform
 fProc->Transform();
 
 Double_t rN=fN;
 
 // Copy the resulting data in the output arrays
 fReOut.Set(fN);
 fImOut.Set(0);
 Double_t re=0;
 for (Int_t i=0; i<fN; i++)
 {
  re=fProc->GetPointReal(i,kFALSE);
  re/=sqrt(rN);
  fReOut.SetAt(re,i);
 }
 
 if (!hist) return;
 
 if ((sel.Contains("Hz") || sel.Contains("t")) && fSample<=0) return;
 
 // Initialize the histogram title
 TString title="";
 if (mode<0) title+="Inverse ";
 title+="DHT ";
 
 TString title2="";
 
 // Define and fill the requested result histogram
 if (sel.Contains("k"))
 {
  hist->SetBins(n,0,n-1);
  title+="index frequency domain";
  if (mode<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Index k;Q[k]";
 }
 if (sel.Contains("f"))
 {
  hist->SetBins(n,0,maxfrac);
  title+="fractional frequency domain";
  if (mode<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Fraction f of sampling rate;Q[f]";
 }
 if (sel.Contains("Hz"))
 {
  hist->SetBins(n,0,maxfrac*fSample);
  title+="actual frequency domain";
  if (mode<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Frequency #nu (Hz);Q[#nu]";
 }
 
 if (sel.Contains("n"))
 {
  hist->SetBins(fN,0,fN);
  title+="sampling time domain";
  if (mode>0) title+=" input";
  if (fSample>0) title2.Form(" (%-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Sample number n;Value X[n]";
 }
 if (sel.Contains("t"))
 {
  hist->SetBins(fN,0,rN/fSample);
  title+="actual time domain";
  if (mode>0) title+=" input";
  if (fSample>0) title2.Form(" (%-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Time t (seconds);Value X[t]";
 }
 
 hist->SetTitle(title);
 
 for (Int_t i=0; i<n; i++)
 {
  if (mode>0) // Forward transform
  {
   if (sel.Contains("n") || sel.Contains("t"))
   {
    hist->SetBinContent(i+1,fReIn.At(i));
   }
   else
   {
    hist->SetBinContent(i+1,fReOut.At(i));
   }
  }
  else // Backward transform
  {
   if (sel.Contains("n") || sel.Contains("t"))
   {
    hist->SetBinContent(i+1,fReOut.At(i));
   }
   else
   {
    hist->SetBinContent(i+1,fReIn.At(i));
   }
  }
 }
}

void NcTransform::Cosine(Int_t type,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 
 if (abs(type)<1 || abs(type)>4 || fN<1 || (abs(type)==1 && fN<2)) return;
 
 // Convert negative type specifications to the corresponding "normal" ones
 Int_t type2=type;
 if (type==-1 || type==-4) type2=abs(type);
 if (type==-2) type2=3;
 if (type==-3) type2=2;
 
 Int_t n=1+fN/2;
 Float_t maxfrac=0.5;
 if (sel.Contains("n") || sel.Contains("t") || sel.Contains("2"))
 {
  n=fN;
  maxfrac=1;
 }
 
 // Comply with the TVirtualFFT conventions
 Int_t kind=type2-1;
 
 // Construct the Fast Fourier Transform (FFT) processor
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 
 fProc=TVirtualFFT::SineCosine(1,&fN,&kind,"ES");
 
 // Enter the input data
 Double_t* data=fReIn.GetArray();
 fProc->SetPoints(data);
 
 // Perform the Discrete Cosine Transform
 fProc->Transform();
 
 Double_t rN=fN;
 
 // Copy the resulting data in the output arrays
 fReOut.Set(fN);
 fImOut.Set(0);
 Double_t re=0;
 for (Int_t i=0; i<fN; i++)
 {
  re=fProc->GetPointReal(i,kFALSE);
  if (type2==1)
  {
   re/=sqrt(2.*(rN-1.));
  }
  else
  {
   re/=sqrt(2.*rN);
  }
  fReOut.SetAt(re,i);
 }
 
 if (!hist) return;
 
 if ((sel.Contains("Hz") || sel.Contains("t")) && fSample<=0) return;
 
 // Initialize the histogram title
 TString title="";
 if (type<0) title+="Inverse ";
 title+="DCT-";
 if (abs(type)==1) title+="I";
 if (abs(type)==2) title+="II";
 if (abs(type)==3) title+="III";
 if (abs(type)==4) title+="IV";
 title+=" ";
 
 TString title2="";
 
 // Define and fill the requested result histogram
 if (sel.Contains("k"))
 {
  hist->SetBins(n,0,n-1);
  title+="index frequency domain";
  if (type<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Index k;Q[k]";
 }
 if (sel.Contains("f"))
 {
  hist->SetBins(n,0,maxfrac);
  title+="fractional frequency domain";
  if (type<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Fraction f of sampling rate;Q[f]";
 }
 if (sel.Contains("Hz"))
 {
  hist->SetBins(n,0,maxfrac*fSample);
  title+="actual frequency domain";
  if (type<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Frequency #nu (Hz);Q[#nu]";
 }
 
 if (sel.Contains("n"))
 {
  hist->SetBins(fN,0,fN);
  title+="sampling time domain";
  if (type>0) title+=" input";
  if (fSample>0) title2.Form(" (%-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Sample number n;Value X[n]";
 }
 if (sel.Contains("t"))
 {
  hist->SetBins(fN,0,rN/fSample);
  title+="actual time domain";
  if (type>0) title+=" input";
  if (fSample>0) title2.Form(" (%-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Time t (seconds);Value X[t]";
 }
 
 hist->SetTitle(title);
 
 // Determine stepsize in fractional sampling frequency
 Double_t fstep=1./(2.*rN);
 if (type==1) fstep=1./(2.*double(fN-1));
 
 Double_t x=0;
 for (Int_t i=0; i<n; i++)
 {
  x=i;
  if (type2==3 || type2==4) x+=0.5;
  x*=fstep;
 
  if (sel.Contains("n") || sel.Contains("t"))
  {
   if (type>0)
   {
    hist->SetBinContent(i+1,fReIn.At(i));
   }
   else
   {
    hist->SetBinContent(i+1,fReOut.At(i));
   }
  }
  else if (sel.Contains("k"))
  {
   if (type>0)
   {
    hist->SetBinContent(i+1,fReOut.At(i));
   }
   else
   {
    hist->SetBinContent(i+1,fReIn.At(i));
   }
  }
  else if (sel.Contains("f"))
  {
   if (type>0)
   {
    hist->Fill(x,fReOut.At(i));
   }
   else
   {
    hist->Fill(x,fReIn.At(i));
   }
  }
  else
  {
   x*=fSample;
   if (type>0)
   {
    hist->Fill(x,fReOut.At(i));
   }
   else
   {
    hist->Fill(x,fReIn.At(i));
   }
  }
 }
}

void NcTransform::Sine(Int_t type,TH1* hist,TString sel)
{
 
 fReOut.Set(0);
 fImOut.Set(0);
 
 if (abs(type)<1 || abs(type)>4 || fN<1 || (abs(type)==1 && fN<2)) return;
 
 // Convert negative type specifications to the corresponding "normal" ones
 Int_t type2=type;
 if (type==-1 || type==-4) type2=abs(type);
 if (type==-2) type2=3;
 if (type==-3) type2=2;
 
 Int_t n=1+fN/2;
 Float_t maxfrac=0.5;
 if (sel.Contains("n") || sel.Contains("t") || sel.Contains("2"))
 {
  n=fN;
  maxfrac=1;
 }
 
 // Comply with the TVirtualFFT conventions
 Int_t kind=type2+3;
 
 // Construct the Fast Fourier Transform (FFT) processor
 if (fProc)
 {
  delete fProc;
  fProc=0;
 }
 
 fProc=TVirtualFFT::SineCosine(1,&fN,&kind,"ES K");
 
 // Enter the input data
 Double_t* data=fReIn.GetArray();
 fProc->SetPoints(data);
 
 // Perform the Discrete Cosine Transform
 fProc->Transform();
 
 Double_t rN=fN;
 
 // Copy the resulting data in the output arrays
 fReOut.Set(fN);
 fImOut.Set(0);
 Double_t re=0;
 for (Int_t i=0; i<fN; i++)
 {
  re=fProc->GetPointReal(i,kFALSE);
  if (type2==1)
  {
   re/=sqrt(2.*(rN+1.));
  }
  else
  {
   re/=sqrt(2.*rN);
  }
  fReOut.SetAt(re,i);
 }
 
 if (!hist) return;
 
 if ((sel.Contains("Hz") || sel.Contains("t")) && fSample<=0) return;
 
 // Initialize the histogram title
 TString title="";
 if (type<0) title+="Inverse ";
 title+="DST-";
 if (abs(type)==1) title+="I";
 if (abs(type)==2) title+="II";
 if (abs(type)==3) title+="III";
 if (abs(type)==4) title+="IV";
 title+=" ";
 
 TString title2="";
 
 // Define and fill the requested result histogram
 if (sel.Contains("k"))
 {
  hist->SetBins(n,0,n-1);
  title+="index frequency domain";
  if (type<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Index k;Q[k]";
 }
 if (sel.Contains("f"))
 {
  hist->SetBins(n,0,maxfrac);
  title+="fractional frequency domain";
  if (type<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Fraction f of sampling rate;Q[f]";
 }
 if (sel.Contains("Hz"))
 {
  hist->SetBins(n,0,maxfrac*fSample);
  title+="actual frequency domain";
  if (type<0)
  {
   title+=" (input)";
  }
  else if (fSample>0)
  {
   title2.Form(" (DAQ: %-.3g samples/sec)",fSample);
   title+=title2;
  }
  title+=";Frequency #nu (Hz);Q[#nu]";
 }
 
 if (sel.Contains("n"))
 {
  hist->SetBins(fN,0,fN);
  title+="sampling time domain";
  if (type>0) title+=" input";
  if (fSample>0) title2.Form(" (%-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Sample number n;Value X[n]";
 }
 if (sel.Contains("t"))
 {
  hist->SetBins(fN,0,rN/fSample);
  title+="actual time domain";
  if (type>0) title+=" input";
  if (fSample>0) title2.Form(" (%-.3g samples/sec)",fSample);
  title+=title2;
  title+=";Time t (seconds);Value X[t]";
 }
 
 hist->SetTitle(title);
 
 // Determine stepsize in fractional sampling frequency
 Double_t fstep=1./(2.*rN);
 if (type==1) fstep=1./(2.*double(fN+1));
 
 Double_t x=0;
 for (Int_t i=0; i<n; i++)
 {
  x=i+1;
  if (type2==3 || type2==4) x-=0.5;
  x*=fstep;
 
  if (sel.Contains("n") || sel.Contains("t"))
  {
   if (type>0)
   {
    hist->SetBinContent(i+1,fReIn.At(i));
   }
   else
   {
    hist->SetBinContent(i+1,fReOut.At(i));
   }
  }
  else if (sel.Contains("k"))
  {
   if (type>0)
   {
    hist->SetBinContent(i+1,fReOut.At(i));
   }
   else
   {
    hist->SetBinContent(i+1,fReIn.At(i));
   }
  }
  else if (sel.Contains("f"))
  {
   if (type>0)
   {
    hist->Fill(x,fReOut.At(i));
   }
   else
   {
    hist->Fill(x,fReIn.At(i));
   }
  }
  else
  {
   x*=fSample;
   if (type>0)
   {
    hist->Fill(x,fReOut.At(i));
   }
   else
   {
    hist->Fill(x,fReIn.At(i));
   }
  }
 }
}

TObject* NcTransform::Clone(const char* name) const
{
 
 NcTransform* q=new NcTransform(*this);
 if (name)
 {
  if (strlen(name)) q->SetName(name);
 }
 return q;
}