From e51acbe0fb83c32912ec4e18d6a9f502bcb03cba Mon Sep 17 00:00:00 2001
From: nemu <nemu>
Date: Thu, 4 Sep 2008 07:13:04 +0000
Subject: [PATCH] added apodization

---
 src/tests/fourier/PMusrFourier.cpp | 116 +++++++++++++++++++++++++++--
 src/tests/fourier/PMusrFourier.h   |  17 +++--
 src/tests/fourier/fourier.cpp      |  31 +++++++-
 3 files changed, 148 insertions(+), 16 deletions(-)

diff --git a/src/tests/fourier/PMusrFourier.cpp b/src/tests/fourier/PMusrFourier.cpp
index 05893629..8b18cea9 100644
--- a/src/tests/fourier/PMusrFourier.cpp
+++ b/src/tests/fourier/PMusrFourier.cpp
@@ -75,6 +75,9 @@ PMusrFourier::PMusrFourier(int dataType, PDoubleVector &data, double timeResolut
   fIn  = 0;
   fOut = 0;
 
+  fApodization = F_APODIZATION_NONE;
+  fFilter      = F_FILTER_NONE;
+
   // if endTime == 0 set is to the last time slot
   if (fEndTime == 0.0)
     fEndTime = (fData.size()-1)*fTimeResolution;
@@ -106,7 +109,7 @@ PMusrFourier::PMusrFourier(int dataType, PDoubleVector &data, double timeResolut
     }
   }
 
-cout << endl << "dB = " << 1.0/(F_GAMMA_BAR_MUON * (fEndTime-fStartTime)) << " (G), Bmax = " << 1.0/(F_GAMMA_BAR_MUON*fTimeResolution) << " (G)" << endl;
+cout << endl << "dB = " << 1.0/(2.0 * F_GAMMA_BAR_MUON * (fEndTime-fStartTime)) << " (G), Bmax = " << 1.0/(2.0 * F_GAMMA_BAR_MUON * fTimeResolution) << " (G)" << endl;
 
 
   // try to estimate N0 and Bkg just out of the raw data
@@ -175,17 +178,18 @@ cout << endl << "in ~PMusrFourier() ..." << endl;
 /**
  * <p>
  *
- * \param apodization
- * \param filter
- * \param estimateN0AndBkg
+ * \param apodizationTag 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod., 4=user apod.
+ * \param filterTag 0=no filter, 1=low pass, 2=band pass, 3=high pass, 4=..., n=user filter
  */
-void PMusrFourier::Transform(int apodization, int filter)
+void PMusrFourier::Transform(unsigned int apodizationTag, unsigned int filterTag)
 {
   if (!fValid)
     return;
 
   if (fDataType == F_SINGLE_HISTO) {
-    PrepareSingleHistoFFTwInputData();
+    PrepareSingleHistoFFTwInputData(apodizationTag, filterTag);
+  } else {
+    PrepareFFTwInputData(apodizationTag, filterTag);
   }
 
 // for test only
@@ -212,7 +216,7 @@ for (int p=-180; p<180; p++) {
   }
   fftw_execute(fFFTwPlan);
 
-  if (p==8) {
+  if (p == 7) {
     for (unsigned int j=0; j<fNoOfBins/2; j++) {
       re.SetBinContent(j+1, fOut[j][0]);
       im.SetBinContent(j+1, fOut[j][1]);
@@ -224,6 +228,9 @@ for (int p=-180; p<180; p++) {
   for (unsigned int i=0; i<fNoOfBins/2; i++) {
     sum += fOut[i][1];
   }
+  if (p == 7) {
+    cout << endl << ">> sum = " << sum << endl;
+  }
   sumHist.SetBinContent(p+181, fabs(sum));
 }
 TFile f("test_out.root", "RECREATE");
@@ -406,7 +413,7 @@ cout << endl << ">> N0/per bin=" << A/(PMUON_LIFETIME*1000.0)*fTimeResolution <<
  * <p>
  *
  */
-void PMusrFourier::PrepareSingleHistoFFTwInputData()
+void PMusrFourier::PrepareSingleHistoFFTwInputData(unsigned int apodizationTag, unsigned int filterTag)
 {
   // 1st fill fIn
   unsigned int start = (unsigned int)(fStartTime/fTimeResolution);
@@ -431,6 +438,9 @@ void PMusrFourier::PrepareSingleHistoFFTwInputData()
   for (unsigned int i=0; i<fNoOfData; i++)
     fIn[i][0] -= fN0;
 
+  // 5th apodize data (if wished)
+  ApodizeData(apodizationTag);
+
 // has to be removed after testing
 TH1F test_in("test_in", "test_in", fNoOfBins,
              fStartTime - fTimeResolution/2.0, fNoOfBins*fTimeResolution + fTimeResolution/2.0);
@@ -440,3 +450,93 @@ TFile f("test_in.root", "RECREATE");
 test_in.Write();
 f.Close();
 }
+
+//--------------------------------------------------------------------------
+// PrepareFFTwInputData
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ */
+void PMusrFourier::PrepareFFTwInputData(unsigned int apodizationTag, unsigned int filterTag)
+{
+  // 1st fill fIn
+  unsigned int start = (unsigned int)(fStartTime/fTimeResolution);
+  for (unsigned int i=0; i<fNoOfData; i++) {
+    fIn[i][0] = fDataRebinned[i+start];
+    fIn[i][1] = 0.0;
+  }
+  for (unsigned int i=fNoOfData; i<fNoOfBins; i++) {
+    fIn[i][0] = 0.0;
+    fIn[i][1] = 0.0;
+  }
+
+  // 2nd apodize data (if wished)
+  ApodizeData(apodizationTag);
+
+// has to be removed after testing
+TH1F test_in("test_in", "test_in", fNoOfBins,
+             fStartTime - fTimeResolution/2.0, fNoOfBins*fTimeResolution + fTimeResolution/2.0);
+for (unsigned int i=0; i<fNoOfBins; i++)
+  test_in.SetBinContent(i, fIn[i][0]);
+TFile f("test_in.root", "RECREATE");
+test_in.Write();
+f.Close();
+}
+
+//--------------------------------------------------------------------------
+// ApodizeData
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ * \param apodizationTag
+ */
+void PMusrFourier::ApodizeData(int apodizationTag) {
+
+  const double cweak[3]   = { 0.384093, -0.087577, 0.703484 };
+  const double cmedium[3] = { 0.152442, -0.136176, 0.983734 };
+  const double cstrong[3] = { 0.045335,  0.554883, 0.399782 };
+
+  double c[5];
+  for (unsigned int i=0; i<5; i++) {
+    c[i] = 0.0;
+  }
+
+  switch (apodizationTag) {
+    case F_APODIZATION_NONE:
+      break;
+    case F_APODIZATION_WEAK:
+      c[0] = cweak[0]+cweak[1]+cweak[2];
+      c[1] = -(cweak[1]+2.0*cweak[2]);
+      c[2] = cweak[2];
+      break;
+    case F_APODIZATION_MEDIUM:
+      c[0] = cmedium[0]+cmedium[1]+cmedium[2];
+      c[1] = -(cmedium[1]+2.0*cmedium[2]);
+      c[2] = cmedium[2];
+      break;
+    case F_APODIZATION_STRONG:
+      c[0] = cstrong[0]+cstrong[1]+cstrong[2];
+      c[1] = -2.0*(cstrong[1]+2.0*cstrong[2]);
+      c[2] = cstrong[1]+6.0*cstrong[2];
+      c[3] = -4.0*cstrong[2];
+      c[4] = cstrong[2];
+      break;
+    case F_APODIZATION_USER:
+      cout << endl << ">> User Apodization not yet implemented, sorry ..." << endl;
+      break;
+    default:
+      cout << endl << ">> **ERROR** User Apodization tag " << apodizationTag << " unknown, sorry ..." << endl;
+      break;
+  }
+
+  double q;
+  for (unsigned int i=0; i<fNoOfData; i++) {
+    q = c[0];
+    for (unsigned int j=1; j<5; j++) {
+      q += c[j] * pow((double)i/(double)fNoOfData, 2.0*(double)j);
+    }
+    fIn[i][0] *= q;
+  }
+}
diff --git a/src/tests/fourier/PMusrFourier.h b/src/tests/fourier/PMusrFourier.h
index c82ba727..8b33f062 100644
--- a/src/tests/fourier/PMusrFourier.h
+++ b/src/tests/fourier/PMusrFourier.h
@@ -39,8 +39,9 @@ using namespace std;
 
 #define F_ESTIMATE_N0_AND_BKG true
 
-#define F_SINGLE_HISTO 0
-#define F_ASYMMETRY    1
+#define F_SINGLE_HISTO_RAW 0
+#define F_SINGLE_HISTO     1
+#define F_ASYMMETRY        2
 
 #define F_APODIZATION_NONE   0
 #define F_APODIZATION_WEAK   1
@@ -87,8 +88,9 @@ class PMusrFourier
 
     virtual void SetN0(double n0) { fN0 = n0; }
     virtual void SetBkg(double bkg) { fBkg = bkg; }
+    virtual void SetPhase(double phase) { fPhase = phase; }
 
-    virtual void Transform(int apodization = -1, int filter = -1);
+    virtual void Transform(unsigned int apodizationTag = 0, unsigned int filter = 0);
 
     virtual double GetFieldResolution() { return fFieldResolution; }
     virtual double GetPhaseCorrection() { return fPhaseCorrection; }
@@ -103,11 +105,12 @@ class PMusrFourier
     bool fValid;
 
     int fDataType; ///< 0=Single Histo, 1=Asymmetry
-    int fUseApodization; ///< 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod., 4=user apod.
-    int fUseFilter; ///< 0=no filter, 1=low pass, 2=band pass, 3=high pass, 4=..., n=user filter
+    int fApodization; ///< 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod., 4=user apod.
+    int fFilter; ///< 0=no filter, 1=low pass, 2=band pass, 3=high pass, 4=..., n=user filter
 
     double fN0;
     double fBkg;
+    double fPhase;
 
     double fTimeResolution;
     double fStartTime;
@@ -126,7 +129,9 @@ class PMusrFourier
     fftw_complex *fIn;
     fftw_complex *fOut;
 
-    virtual void PrepareSingleHistoFFTwInputData();
+    virtual void PrepareSingleHistoFFTwInputData(unsigned int apodizationTag, unsigned int filterTag);
+    virtual void PrepareFFTwInputData(unsigned int apodizationTag, unsigned int filterTag);
+    virtual void ApodizeData(int apodizationTag);
     virtual void EstimateN0AndBkg();
 };
 
diff --git a/src/tests/fourier/fourier.cpp b/src/tests/fourier/fourier.cpp
index ecc28b23..33972669 100644
--- a/src/tests/fourier/fourier.cpp
+++ b/src/tests/fourier/fourier.cpp
@@ -159,11 +159,38 @@ cout << endl << "#bins=" << histo->GetNbinsX();
     cin >> zeroPaddingPower;
   }
 
+  bool estimateN0AndBkg = true;
+  double N0, bkg, phase;
+  cout << endl << ">> Do you want to give N0, Bkg, and phase explicitly (y/n)? ";
+  cin >> answer;
+  if (strstr(answer, "y")) {
+    estimateN0AndBkg = false;
+    cout << endl << ">> N0    = ";
+    cin >> N0;
+    cout << endl << ">> bkg   = ";
+    cin >> bkg;
+    cout << endl << ">> phase = ";
+    cin >> phase;
+  }
+
   PMusrFourier fourier(F_SINGLE_HISTO, data, timeResolution, startTime, endTime, 
-                       rebin, zeroPaddingPower, F_ESTIMATE_N0_AND_BKG);
+                       rebin, zeroPaddingPower, estimateN0AndBkg);
 
   if (fourier.IsValid()) {
-    fourier.Transform(F_APODIZATION_NONE, F_FILTER_NONE);
+    if (!estimateN0AndBkg) {
+      fourier.SetN0(N0);
+      fourier.SetBkg(bkg);
+      fourier.SetPhase(phase);
+    }
+
+    cout << endl << ">> Do you wish to apodize your data (y/n)? ";
+    cin >> answer;
+    unsigned int apodizationTag=0;
+    if (strstr(answer, "y")) {
+      cout << endl << ">> apodization (1=weak, 2=medium, 3=strong, 4=user) = ";
+      cin >> apodizationTag;
+    }
+    fourier.Transform(apodizationTag, F_FILTER_NONE);
   } else {
     cout << endl << "**ERROR** something went wrong when invoking the PMusrFourier object ...";
     cout << endl << endl;