Index: trunk/MagicSoft/Mars/mtools/MFFT.cc
===================================================================
--- trunk/MagicSoft/Mars/mtools/MFFT.cc	(revision 3053)
+++ trunk/MagicSoft/Mars/mtools/MFFT.cc	(revision 3054)
@@ -94,5 +94,5 @@
 }
 
-void MFFT::TransformF(const Int_t isign)
+void MFFT::TransformF(const Int_t isign, TArrayF &data)
 {
   
@@ -111,6 +111,6 @@
   for (i=1;i<n;i+=2) {
     if (j > i) {
-      Swap(fDataF[j-1],fDataF[i-1]);
-      Swap(fDataF[j],fDataF[i]);
+      Swap(data[j-1],data[i-1]);
+      Swap(data[j],data[i]);
     }
     m=nn;
@@ -146,10 +146,10 @@
         //
         j          = i+mmax;
-        tempr      = wr*fDataF[j-1] - wi*fDataF[j];
-        tempi      = wr*fDataF[j]   + wi*fDataF[j-1];
-        fDataF[j-1] = fDataF[i-1]   - tempr;
-        fDataF[j]   = fDataF[i]     - tempi;
-        fDataF[i-1] += tempr;
-        fDataF[i]   += tempi;
+        tempr      = wr*data[j-1] - wi*data[j];
+        tempi      = wr*data[j]   + wi*data[j-1];
+        data[j-1] = data[i-1]   - tempr;
+        data[j]   = data[i]     - tempi;
+        data[i-1] += tempr;
+        data[i]   += tempi;
       }
 
@@ -166,5 +166,5 @@
 
 
-void MFFT::TransformD(const Int_t isign)
+void MFFT::TransformD(const Int_t isign, TArrayD &data)
 {
   
@@ -183,6 +183,6 @@
   for (i=1;i<n;i+=2) {
     if (j > i) {
-      Swap(fDataD[j-1],fDataD[i-1]);
-      Swap(fDataD[j],fDataD[i]);
+      Swap(data[j-1],data[i-1]);
+      Swap(data[j],data[i]);
     }
     m=nn;
@@ -218,10 +218,10 @@
         //
         j          = i+mmax;
-        tempr      = wr*fDataD[j-1] - wi*fDataD[j];
-        tempi      = wr*fDataD[j]   + wi*fDataD[j-1];
-        fDataD[j-1] = fDataD[i-1]   - tempr;
-        fDataD[j]   = fDataD[i]     - tempi;
-        fDataD[i-1] += tempr;
-        fDataD[i]   += tempi;
+        tempr      = wr*data[j-1] - wi*data[j];
+        tempi      = wr*data[j]   + wi*data[j-1];
+        data[j-1] = data[i-1]   - tempr;
+        data[j]   = data[i]     - tempi;
+        data[i-1] += tempr;
+        data[i]   += tempi;
       }
 
@@ -262,5 +262,5 @@
     {
       c2    = -0.5;
-      TransformF(1);
+      TransformF(1,fDataF);
     }
   else         // set up backward transform
@@ -325,5 +325,5 @@
       // The inverse transform for the case isign = -1
       //
-      TransformF(-1);  
+      TransformF(-1,fDataF);  
 
       //
@@ -349,5 +349,5 @@
     {
       c2    = -0.5;
-      TransformD(1);
+      TransformD(1,fDataD);
     }
   else         // set up backward transform
@@ -412,5 +412,5 @@
       // The inverse transform for the case isign = -1
       //
-      TransformD(-1);  
+      TransformD(-1,fDataD);  
       
       //
@@ -648,9 +648,109 @@
 }
 
-
-//
-// Power Spectrum Density calculation
-//
-TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist)
+//
+// Power Spectrum Density calculation for TArrayF
+//
+TArrayF* MFFT::PowerSpectrumDensity(const TArrayF *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] = 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=0;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)
+{
+  
+  fDim = array->GetSize();
+  CheckDim(fDim);
+
+  fDataD.Set(fDim);
+  //
+  // Copy the hist into an array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataD[i] = array->At(i);
+
+  RealFTD(1);
+
+  const Int_t dim2  = fDim*fDim;
+  const Int_t dim05 = fDim/2;
+  Float_t c02;
+  Float_t ck2;
+  Float_t cn2;
+  
+  TArrayD *newarray = new TArrayD(dim05);
+
+  //
+  // Fill the new histogram: 
+  //
+  // 1) P(0) = 1/(N*N) |C(0)|*|C(0)|
+  //
+  c02 = (fDataD[0]*fDataD[0]);
+  //  newarray->AddAt(c02/dim2,0);
+  //
+  // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
+  //
+  for (Int_t k=0;k<dim05-1;k++)
+    {
+      const Int_t k2 = k+k;
+      ck2 = (fDataD[k2]*fDataD[k2] + fDataD[k2+1]*fDataD[k2+1]);
+      newarray->AddAt(ck2/dim2,k);
+    }
+  //
+  // 3) P(N) = 1/(N*N) (|C(n/2)|*|C(n/2)|)
+  //
+  cn2 = (fDataD[1]*fDataD[1]);
+  //  newarray->AddAt(cn2,dim05-1);
+  
+  return newarray;
+}
+
+
+//
+// Power Spectrum Density calculation for TH1
+//
+TH1F* MFFT::PowerSpectrumDensity(const TH1 *hist)
 {
 
@@ -676,5 +776,5 @@
   //
   c02 = (fDataF[0]*fDataF[0]);
-  newhist->Fill(c02/dim2);
+  newhist->Fill(0.,c02/dim2);
   //
   // 2) P(k) = 1/(N*N) (|C(k)|*|C(k)|))
@@ -682,7 +782,6 @@
   for (Int_t k=2;k<fDim;k+=2)
     {
-      Int_t ki  = k+1;
-      ck2 = (fDataF[k]*fDataF[k] + fDataF[ki]*fDataF[ki]);
-      newhist->Fill(ck2/dim2);
+      ck2 = (fDataF[k]*fDataF[k] + fDataF[k+1]*fDataF[k+1]);
+      newhist->Fill(k/2.,ck2/dim2);
     }
   //
@@ -690,7 +789,24 @@
   //
   cn2 = (fDataF[1]*fDataF[1]);
-  newhist->Fill(cn2/dim2);
+  newhist->Fill(fDim/2.-1.,cn2/dim2);
   
   return newhist;
+}
+
+
+//
+// Power Spectrum Density calculation for TH1I
+//
+TH1F* MFFT::PowerSpectrumDensity(const TH1F *hist)
+{
+  return PowerSpectrumDensity((TH1*)hist);
+}
+
+//
+// Power Spectrum Density calculation for TH1I
+//
+TH1F* MFFT::PowerSpectrumDensity(const TH1I *hist)
+{
+  return PowerSpectrumDensity((TH1*)hist);
 }
 
@@ -721,9 +837,4 @@
 {
   
-  TString name = hist->GetName();
-  name += " PSD";
-  TString title = hist->GetTitle();
-  title += " - Power Spectrum Density";
-
   // number of entries
   fDim = hist->GetNbinsX();
@@ -735,5 +846,5 @@
   // Nyquist frequency
   Axis_t fcrit = 1./(2.*delta);
-  Axis_t low = 0.;
+  Axis_t low = -0.5;
   Axis_t up  = fcrit;
 
@@ -741,12 +852,81 @@
     {
     case 0: 
-      return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
+      return new TH1F(Form("%s%s",hist->GetName()," PSD"),
+                      Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
+                      fDim/2,low,up);
       break;
     case 1:
-      return new TH1D(name.Data(),title.Data(),fDim/2+1,low,up);
+      return new TH1D(Form("%s%s",hist->GetName()," PSD"),
+                      Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
+                      fDim/2,low,up);
       break;
     default:
-      return new TH1F(name.Data(),title.Data(),fDim/2+1,low,up);
+      return new TH1F(Form("%s%s",hist->GetName()," PSD"),
+                      Form("%s%s",hist->GetTitle()," - Power Spectrum Density"),
+                      fDim/2,low,up);
       break;
     }
 }
+
+//
+// Real function spectrum with data windowing
+//
+TArrayF* MFFT::RealFunctionSpectrum(const TArrayF *data) 
+{
+  
+  fDim = data->GetSize();
+  CheckDim(fDim);
+
+  fDataF.Set(fDim);
+  //
+  // Copy the hist into an array
+  //
+  for (Int_t i=0;i<fDim;i++)
+    fDataF[i] = data->At(i);
+
+  fWindowF.Set(fDim);
+
+  Int_t dim2 = fDim/2;
+
+  TArrayF *power = new TArrayF(dim2);
+
+  // 
+  // Start program spctrm from NR's
+  //
+  Float_t w, facp, facm, sumw=0.;
+  
+  facm = dim2;
+  facp = 1./dim2;
+  
+  for (Int_t j=0;j<dim2;j++)
+    {
+      Int_t j2 = j+j;
+      w     = ApplyWindow(j,facm,facp);
+      sumw += w*w;
+      fWindowF[j2]   = fDataF[j]*w;
+      fWindowF[j2+1] = fDataF[dim2+j]*w;
+    }
+  
+  TransformF(1,fWindowF);
+  
+  power->AddAt(fWindowF[0]*fWindowF[0] + fWindowF[1]*fWindowF[1],0);
+
+  //  power->AddAt(fWindowF[0]*fWindowF[0],0);
+  //  power->AddAt(fWindowF[1]*fWindowF[1],dim2-1);  
+
+
+  for (Int_t j=1;j<dim2;j++)
+    //  for (Int_t j=1;j<dim2;j++)
+    {
+      Int_t j2 = j+j;
+      Float_t buf = fWindowF[j2+1]     *fWindowF[j2+1] 
+                  + fWindowF[j2  ]     *fWindowF[j2  ] 
+                  + fWindowF[fDim-j2+1]*fWindowF[fDim-j2+1] 
+                  + fWindowF[fDim-j2  ]*fWindowF[fDim-j2  ] ;
+      power->AddAt(buf/sumw/(fDim+fDim),j);
+    }
+  
+  return power;
+
+}
+
Index: trunk/MagicSoft/Mars/mtools/MSimulatedAnnealing.cc
===================================================================
--- trunk/MagicSoft/Mars/mtools/MSimulatedAnnealing.cc	(revision 3053)
+++ trunk/MagicSoft/Mars/mtools/MSimulatedAnnealing.cc	(revision 3054)
@@ -45,5 +45,5 @@
 //  2) a TArrayF y(ndim+1)                                                  //
 //     whose components must be pre-initialized to the values of            //
-//      FunctionToMinimize evaluated at the fNdim+1 vertices (rows) of fP   //
+//      FunctionToMinimize evaluated at the fNdim+1 vertices (rows) of p    //
 //  3) a TArrayF p0(ndim)                                                   //
 //     whose components contain the lower simplex borders                   //