Changeset 3125 for trunk/MagicSoft/Mars/mtools/MFFT.cc

Timestamp:

02/12/04 19:56:16 (21 years ago)

Author:

gaug

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mtools/MFFT.cc (modified) (1 diff)

trunk/MagicSoft/Mars/mtools/MFFT.cc

@@ -700 +700 @@
 
 
+//
+// Power Spectrum Density calculation for TArrayI
+//
+TArrayF* MFFT::PowerSpectrumDensity(const TArrayI *array)
+{
+
+  fDim = array->GetSize();
+  CheckDim(fDim);
+
+  fDataF.Set(fDim);
+  //
+  // Copy the hist into an array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataF[i] = (Float_t)array->At(i);
+
+  RealFTF(1);
+
+  const Int_t dim2  = fDim*fDim;
+  const Int_t dim05 = fDim/2;
+  Float_t c02;
+  Float_t ck2;
+  Float_t cn2;
+  
+  TArrayF *newarray = new TArrayF(dim05);
+
+  //
+  // Fill the new histogram: 
+  //
+  // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
+  //
+  c02 = (fDataF[0]*fDataF[0]);
+  //  newarray->AddAt(c02/dim2,0);
+  //
+  // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
+  //
+  for (Int_t k=1;k<dim05-1;k++)
+    {
+      const Int_t k2 = k+k;
+      ck2 = (fDataF[k2]*fDataF[k2] + fDataF[k2+1]*fDataF[k2+1]);
+      newarray->AddAt(ck2/dim2,k);
+    }
+  //
+  // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
+  //
+  cn2 = (fDataF[1]*fDataF[1]);
+  //  newarray->AddAt(cn2,dim05-1);
+  
+  return newarray;
+}
+
+
 TArrayD* MFFT::PowerSpectrumDensity(const TArrayD *array)
 {