From b8b7b7665f868a0bd1b86558ffc33af1a286bb10 Mon Sep 17 00:00:00 2001
From: Suter Andreas <andreas.suter@psi.ch>
Date: Fri, 15 Jan 2016 18:26:47 +0100
Subject: [PATCH] some more work towards asymmetry RRF

---
 src/classes/PRunAsymmetryRRF.cpp   | 27 ++++++++++++++++++++-------
 src/classes/PRunListCollection.cpp |  6 +++---
 src/musrfit.cpp                    | 26 ++++++++++++++++++++++++++
 3 files changed, 49 insertions(+), 10 deletions(-)

diff --git a/src/classes/PRunAsymmetryRRF.cpp b/src/classes/PRunAsymmetryRRF.cpp
index bd21a5b8..eef4df0c 100644
--- a/src/classes/PRunAsymmetryRRF.cpp
+++ b/src/classes/PRunAsymmetryRRF.cpp
@@ -386,6 +386,11 @@ void PRunAsymmetryRRF::SetFitRangeBin(const TString fitRange)
  */
 void PRunAsymmetryRRF::CalcNoOfFitBins()
 {
+  cout << "debug> fData.GetValue()->size()=" << fData.GetValue()->size() << endl;
+  cout << "debug> fFitStartTime=" << fFitStartTime << endl;
+  cout << "debug> fData.GetDataTimeStart()=" << fData.GetDataTimeStart() << ", fData.GetDataTimeStep()=" << fData.GetDataTimeStep() << endl;
+  cout << "debug> -----" << endl;
+
   // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly
   Int_t startTimeBin = static_cast<Int_t>(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep()));
   if (startTimeBin < 0)
@@ -867,9 +872,7 @@ Bool_t PRunAsymmetryRRF::PrepareFitData()
   Int_t lgb = fgb + lgb_offset;
   Int_t dt0 = (Int_t)fT0s[0]-(Int_t)fT0s[1];
 
-  cout << "debug> fGoodBins[0]=" << fGoodBins[0] << ", fGoodBins[2]=" << fGoodBins[2] << ", fT0s[0]=" << fT0s[0] << ", fT0s[1]=" << fT0s[1] << endl;
-  cout << "debug> fgbOffset=" << fgbOffset << endl;
-  cout << "debug> fgb=" << fgb << ", lgb=" << lgb << ", dt0=" << dt0 << endl;
+  cout << "debug> fgb=" << fgb << ", lgb=" << lgb << endl;
 
   PDoubleVector asym;
   PDoubleVector asymErr;
@@ -885,7 +888,7 @@ Bool_t PRunAsymmetryRRF::PrepareFitData()
     asym.push_back(asymVal);
     eff = fForwardErr[i];
     ebb = fBackwardErr[i-dt0];
-    if (asymVal != 0.0)
+    if ((asymVal != 0.0) && (ff+bb) > 0.0)
       asymValErr = sqrt(2)/(ff+bb)*sqrt(bb*eff+ff*ebb);
     else
       asymValErr = 1.0;
@@ -897,13 +900,15 @@ Bool_t PRunAsymmetryRRF::PrepareFitData()
   Double_t wRRF = globalBlock->GetRRFFreq("Mc");
   Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0;
 
-  Double_t startTime = fTimeResolution * (fT0s[0]+(Double_t)fgbOffset);
+  Double_t startTime = fTimeResolution * (Double_t)fgbOffset;
   Double_t time=0.0;
   for (UInt_t i=0; i<asym.size(); i++) {
     time = startTime + i*fTimeResolution;
     asym[i] *= 2.0*cos(wRRF*time+phaseRRF);
   }
 
+  cout << "debug> before packing: startTime=" << startTime << ", endTime=" << startTime+asym.size()*fTimeResolution << endl;
+
   // 3rd: rrf packing
   PDoubleVector asymRRF;
   asymVal = 0.0;
@@ -926,11 +931,19 @@ Bool_t PRunAsymmetryRRF::PrepareFitData()
     asymValErr += asymErr[i]*asymErr[i];
   }
 
+  cout << "debug> after packing: startTime=" << startTime + fTimeResolution*((Double_t)(fRRFPacking-1)/2.0) << ", endTime=" << startTime + fTimeResolution*((Double_t)(fRRFPacking-1)/2.0)+asymRRF.size()*fTimeResolution*(Double_t)(fRRFPacking) << endl;
 
-  fData.SetDataTimeStart(startTime+(Double_t)(fRRFPacking-1)/2.0);
+  ofstream fout("_data.dat", ofstream::out);
+  for (UInt_t i=0; i<asymRRF.size(); i++) {
+    fout << startTime + fTimeResolution*((Double_t)(fRRFPacking-1)/2.0) + i*fTimeResolution*(Double_t)fRRFPacking << ", " << asymRRF[i] << ", " << asymRRFErr[i] << endl;
+  }
+  fout.close();
+
+
+  fData.SetDataTimeStart(startTime+fTimeResolution*((Double_t)(fRRFPacking-1)/2.0));
   fData.SetDataTimeStep(fTimeResolution*(Double_t)fRRFPacking);
 
-  for (UInt_t i=0; i<asym.size(); i++) {
+  for (UInt_t i=0; i<asymRRF.size(); i++) {
     fData.AppendValue(asymRRF[i]);
     fData.AppendErrorValue(asymRRFErr[i]);
   }
diff --git a/src/classes/PRunListCollection.cpp b/src/classes/PRunListCollection.cpp
index c5f57396..c9c5e56f 100644
--- a/src/classes/PRunListCollection.cpp
+++ b/src/classes/PRunListCollection.cpp
@@ -832,11 +832,11 @@ PRunData* PRunListCollection::GetAsymmetryRRF(UInt_t index, EDataSwitch tag)
         return 0;
       }
 
-      fRunAsymmetryList[index]->CalcTheory();
-      data = fRunAsymmetryList[index]->GetData();
+      fRunAsymmetryRRFList[index]->CalcTheory();
+      data = fRunAsymmetryRRFList[index]->GetData();
       break;
     case kRunNo: // called from PMusrCanvas
-      for (UInt_t i=0; i<fRunAsymmetryList.size(); i++) {
+      for (UInt_t i=0; i<fRunAsymmetryRRFList.size(); i++) {
         if (fRunAsymmetryRRFList[i]->GetRunNo() == index) {
           data = fRunAsymmetryRRFList[i]->GetData();
           break;
diff --git a/src/musrfit.cpp b/src/musrfit.cpp
index 13e72c18..f5f5ac44 100644
--- a/src/musrfit.cpp
+++ b/src/musrfit.cpp
@@ -216,6 +216,19 @@ void musrfit_dump_ascii(char *fileName, PRunListCollection *runList)
     }
   }
 
+  // asymmetry RRF
+  size = runList->GetNoOfAsymmetryRRF();
+  if (size > 0) {
+    for (unsigned int i=0; i<size; i++) {
+      data = runList->GetAsymmetryRRF(i);
+      if (data) {
+        // dump data
+        musrfit_write_ascii(fln, data, runCounter);
+        runCounter++;
+      }
+    }
+  }
+
   // muMinus
   size = runList->GetNoOfMuMinus();
   if (size > 0) {
@@ -362,6 +375,19 @@ void musrfit_dump_root(char *fileName, PRunListCollection *runList)
     }
   }
 
+  // asymmetry RRF
+  size = runList->GetNoOfAsymmetryRRF();
+  if (size > 0) {
+    for (unsigned int i=0; i<size; i++) {
+      data = runList->GetAsymmetryRRF(i);
+      if (data) {
+        // dump data
+        musrfit_write_root(f, fln, data, runCounter);
+        runCounter++;
+      }
+    }
+  }
+
   // muMinus
   size = runList->GetNoOfMuMinus();
   if (size > 0) {