first implementation of the single histo RRF fit (no plotting yet). Seems to work but the error estimate of the RRF asymmetry seems to be too large. Needs to be checked.

src/classes/Makefile.am

@@ -21,6 +21,7 @@ h_sources = \
   ../include/PRunNonMusr.h \
   ../include/PRunMuMinus.h \
   ../include/PRunSingleHisto.h \
+  ../include/PRunSingleHistoRRF.h \
   ../include/PStartupHandler.h \
   ../include/PTheory.h
 
@@ -65,6 +66,7 @@ cpp_sources = \
   PRunNonMusr.cpp \
   PRunMuMinus.cpp \
   PRunSingleHisto.cpp \
+  PRunSingleHistoRRF.cpp \
   PStartupHandler.cpp \
   PTheory.cpp
 

src/classes/PFitterFcn.cpp

@@ -77,11 +77,13 @@ Double_t PFitterFcn::operator()(const std::vector<Double_t>& par) const
 
   if (fUseChi2) { // chi square
     value += fRunListCollection->GetSingleHistoChisq(par);
+    value += fRunListCollection->GetSingleHistoRRFChisq(par);
     value += fRunListCollection->GetAsymmetryChisq(par);
     value += fRunListCollection->GetMuMinusChisq(par);
     value += fRunListCollection->GetNonMusrChisq(par);
   } else { // max likelihood
     value += fRunListCollection->GetSingleHistoMaximumLikelihood(par);
+    value += fRunListCollection->GetSingleHistoRRFMaximumLikelihood(par);
     value += fRunListCollection->GetAsymmetryMaximumLikelihood(par);
     value += fRunListCollection->GetMuMinusMaximumLikelihood(par);
     value += fRunListCollection->GetNonMusrMaximumLikelihood(par);

src/classes/PMsrHandler.cpp

@@ -636,7 +636,13 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages)
           fout.width(16);
           fout << left << "rrf_freq ";
           fout.width(8);
-          fout << left << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data());
+
+cout << "debug> PMsrHandler::WriteMsrLogFile(): fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data())=" << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) << endl;
+neededPrec = LastSignificant(fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()),10);
+cout << "debug> PMsrHandler::WriteMsrLogFile(): neededPrec=" << neededPrec << endl;
+fout.precision(neededPrec);
+
+          fout << left << fixed << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data());
           fout << " " << fGlobal.GetRRFUnit();
           fout << endl;
         } else if (sstr.BeginsWith("rrf_phase", TString::kIgnoreCase) && (fGlobal.GetFitType() == MSR_FITTYPE_SINGLE_HISTO_RRF)) {
@@ -1650,6 +1656,7 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map<UInt_t, TString> *co
       fout.width(16);
       fout << left << "rrf_freq ";
       fout.width(8);
+cout << "debug> PMsrHandler::WriteMsrFile(): fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data())=" << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data()) << endl;
       fout << left << fGlobal.GetRRFFreq(fGlobal.GetRRFUnit().Data());
       fout << " " << fGlobal.GetRRFUnit();
       fout << endl;

src/classes/PMusr.cpp

@@ -752,21 +752,25 @@ Double_t PMsrGlobalBlock::GetRRFFreq(const char *unit)
     return freq;
   }
 
+cout << endl << "debug> PMsrGlobalBlock::GetRRFFreq(" << unit << "): unitTag=" << unitTag << ", fRRFFreq=" << fRRFFreq;
+
   // calc the conversion factor
   if (unitTag == fRRFUnitTag)
     freq = fRRFFreq;
   else if ((unitTag == RRF_UNIT_MHz) && (fRRFUnitTag == RRF_UNIT_Mcs))
     freq = fRRFFreq/TMath::TwoPi();
   else if ((unitTag == RRF_UNIT_MHz) && (fRRFUnitTag == RRF_UNIT_T))
-    freq = fRRFFreq*GAMMA_BAR_MUON;
+    freq = fRRFFreq*1e4*GAMMA_BAR_MUON; // 1e4 need for T -> G since GAMMA_BAR_MUON is given in MHz/G
   else if ((unitTag == RRF_UNIT_Mcs) && (fRRFUnitTag == RRF_UNIT_MHz))
     freq = fRRFFreq*TMath::TwoPi();
   else if ((unitTag == RRF_UNIT_Mcs) && (fRRFUnitTag == RRF_UNIT_T))
-    freq = fRRFFreq*TMath::TwoPi()*GAMMA_BAR_MUON;
+    freq = fRRFFreq*1e4*TMath::TwoPi()*GAMMA_BAR_MUON; // 1e4 need for T -> G since GAMMA_BAR_MUON is given in MHz/G
   else if ((unitTag == RRF_UNIT_T) && (fRRFUnitTag == RRF_UNIT_MHz))
-    freq = fRRFFreq/GAMMA_BAR_MUON;
+    freq = fRRFFreq/GAMMA_BAR_MUON*1e-4; // 1e-4 need for G -> T since GAMMA_BAR_MUON is given in MHz/G
   else if ((unitTag == RRF_UNIT_T) && (fRRFUnitTag == RRF_UNIT_Mcs))
-    freq = fRRFFreq/(TMath::TwoPi()*GAMMA_BAR_MUON);
+    freq = fRRFFreq/(TMath::TwoPi()*GAMMA_BAR_MUON)*1e-4; // 1e-4 need for G -> T since GAMMA_BAR_MUON is given in MHz/G
+
+cout << endl << "debug> PMsrGlobalBlock::GetRRFFreq(" << unit << "): freq=" << freq << endl;
 
   return freq;
 }

src/classes/PRunListCollection.cpp

@@ -58,6 +58,12 @@ PRunListCollection::~PRunListCollection()
   }
   fRunSingleHistoList.clear();
 
+  for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++) {
+    fRunSingleHistoRRFList[i]->CleanUp();
+    fRunSingleHistoRRFList[i]->~PRunSingleHistoRRF();
+  }
+  fRunSingleHistoRRFList.clear();
+
   for (UInt_t i=0; i<fRunAsymmetryList.size(); i++) {
     fRunAsymmetryList[i]->CleanUp();
     fRunAsymmetryList[i]->~PRunAsymmetry();
@@ -100,12 +106,20 @@ Bool_t PRunListCollection::Add(Int_t runNo, EPMusrHandleTag tag)
     fitType = (*fMsrInfo->GetMsrGlobal()).GetFitType();
   }
 
+cout << "debug> PRunListCollection::Add(): (runNo: " << runNo << "), fitType = " << fitType << endl;
+
   switch (fitType) {
     case PRUN_SINGLE_HISTO:
       fRunSingleHistoList.push_back(new PRunSingleHisto(fMsrInfo, fData, runNo, tag));
       if (!fRunSingleHistoList[fRunSingleHistoList.size()-1]->IsValid())
         success = false;
       break;
+    case PRUN_SINGLE_HISTO_RRF:
+cout << "debug> PRunListCollection::Add(): add RRF single histo run to PRunListCollection (runNo: " << runNo << ")" << endl;
+      fRunSingleHistoRRFList.push_back(new PRunSingleHistoRRF(fMsrInfo, fData, runNo, tag));
+      if (!fRunSingleHistoRRFList[fRunSingleHistoRRFList.size()-1]->IsValid())
+        success = false;
+      break;
     case PRUN_ASYMMETRY:
       fRunAsymmetryList.push_back(new PRunAsymmetry(fMsrInfo, fData, runNo, tag));
       if (!fRunAsymmetryList[fRunAsymmetryList.size()-1]->IsValid())
@@ -147,6 +161,8 @@ void PRunListCollection::SetFitRange(const TString fitRange)
 {
   for (UInt_t i=0; i<fRunSingleHistoList.size(); i++)
     fRunSingleHistoList[i]->SetFitRangeBin(fitRange);
+  for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
+    fRunSingleHistoRRFList[i]->SetFitRangeBin(fitRange);
   for (UInt_t i=0; i<fRunAsymmetryList.size(); i++)
     fRunAsymmetryList[i]->SetFitRangeBin(fitRange);
   for (UInt_t i=0; i<fRunMuMinusList.size(); i++)
@@ -169,6 +185,8 @@ void PRunListCollection::SetFitRange(const PDoublePairVector fitRange)
 {
   for (UInt_t i=0; i<fRunSingleHistoList.size(); i++)
     fRunSingleHistoList[i]->SetFitRange(fitRange);
+  for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
+    fRunSingleHistoRRFList[i]->SetFitRange(fitRange);
   for (UInt_t i=0; i<fRunAsymmetryList.size(); i++)
     fRunAsymmetryList[i]->SetFitRange(fitRange);
   for (UInt_t i=0; i<fRunMuMinusList.size(); i++)
@@ -198,6 +216,27 @@ Double_t PRunListCollection::GetSingleHistoChisq(const std::vector<Double_t>& pa
   return chisq;
 }
 
+//--------------------------------------------------------------------------
+// GetSingleHistoRRFChisq (public)
+//--------------------------------------------------------------------------
+/**
+ * <p>Calculates chi-square of <em>all</em> single histogram RRF runs of a msr-file.
+ *
+ * <b>return:</b>
+ * - chi-square of all single histogram RRF runs of the msr-file
+ *
+ * \param par fit parameter vector
+ */
+Double_t PRunListCollection::GetSingleHistoRRFChisq(const std::vector<Double_t>& par) const
+{
+  Double_t chisq = 0.0;
+
+  for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
+    chisq += fRunSingleHistoRRFList[i]->CalcChiSquare(par);
+
+  return chisq;
+}
+
 //--------------------------------------------------------------------------
 // GetAsymmetryChisq (public)
 //--------------------------------------------------------------------------
@@ -299,6 +338,9 @@ Double_t PRunListCollection::GetSingleHistoChisqExpected(const std::vector<Doubl
   case PRUN_SINGLE_HISTO:
     expectedChisq = fRunSingleHistoList[subIdx]->CalcChiSquareExpected(par);
     break;
+  case PRUN_SINGLE_HISTO_RRF:
+    expectedChisq = fRunSingleHistoRRFList[subIdx]->CalcChiSquareExpected(par);
+    break;
   case PRUN_ASYMMETRY:
     expectedChisq = fRunAsymmetryList[subIdx]->CalcChiSquareExpected(par);
     break;
@@ -353,6 +395,9 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector<Double_t>& par,
   case PRUN_SINGLE_HISTO:
     chisq = fRunSingleHistoList[subIdx]->CalcChiSquare(par);
     break;
+  case PRUN_SINGLE_HISTO_RRF:
+    chisq = fRunSingleHistoRRFList[subIdx]->CalcChiSquare(par);
+    break;
   case PRUN_ASYMMETRY:
     chisq = fRunAsymmetryList[subIdx]->CalcChiSquare(par);
     break;
@@ -376,7 +421,7 @@ Double_t PRunListCollection::GetSingleRunChisq(const std::vector<Double_t>& par,
  * <p>Calculates log max-likelihood of <em>all</em> single histogram runs of a msr-file.
  *
  * <b>return:</b>
- * - chi-square of all single histogram runs of the msr-file
+ * - log max-likelihood of all single histogram runs of the msr-file
  *
  * \param par fit parameter vector
  */
@@ -390,6 +435,27 @@ Double_t PRunListCollection::GetSingleHistoMaximumLikelihood(const std::vector<D
   return mlh;
 }
 
+//--------------------------------------------------------------------------
+// GetSingleHistoRRFMaximumLikelihood (public)
+//--------------------------------------------------------------------------
+/**
+ * <p>Calculates log max-likelihood of <em>all</em> single histogram RRF runs of a msr-file.
+ *
+ * <b>return:</b>
+ * - log max-likelihood of all single histogram runs of the msr-file
+ *
+ * \param par fit parameter vector
+ */
+Double_t PRunListCollection::GetSingleHistoRRFMaximumLikelihood(const std::vector<Double_t>& par) const
+{
+  Double_t mlh = 0.0;
+
+  for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
+    mlh += fRunSingleHistoRRFList[i]->CalcMaxLikelihood(par);
+
+  return mlh;
+}
+
 //--------------------------------------------------------------------------
 // GetAsymmetryMaximumLikelihood (public)
 //--------------------------------------------------------------------------
@@ -419,7 +485,7 @@ Double_t PRunListCollection::GetAsymmetryMaximumLikelihood(const std::vector<Dou
  * <p>Calculates log max-likelihood of <em>all</em> mu minus runs of a msr-file.
  *
  * <b>return:</b>
- * - chi-square of all mu minus runs of the msr-file
+ * - log max-likelihood of all mu minus runs of the msr-file
  *
  * \param par fit parameter vector
  */
@@ -493,6 +559,9 @@ UInt_t PRunListCollection::GetNoOfBinsFitted(const UInt_t idx) const
   case PRUN_SINGLE_HISTO:
     result = fRunSingleHistoList[subIdx]->GetNoOfFitBins();
     break;
+  case PRUN_SINGLE_HISTO_RRF:
+    result = fRunSingleHistoRRFList[subIdx]->GetNoOfFitBins();
+    break;
   case PRUN_ASYMMETRY:
     result = fRunAsymmetryList[subIdx]->GetNoOfFitBins();
     break;
@@ -526,6 +595,9 @@ UInt_t PRunListCollection::GetTotalNoOfBinsFitted() const
   for (UInt_t i=0; i<fRunSingleHistoList.size(); i++)
     counts += fRunSingleHistoList[i]->GetNoOfFitBins();
 
+  for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++)
+    counts += fRunSingleHistoRRFList[i]->GetNoOfFitBins();
+
   for (UInt_t i=0; i<fRunAsymmetryList.size(); i++)
     counts += fRunAsymmetryList[i]->GetNoOfFitBins();
 
@@ -581,6 +653,49 @@ PRunData* PRunListCollection::GetSingleHisto(UInt_t index, EDataSwitch tag)
   return data;
 }
 
+//--------------------------------------------------------------------------
+// GetSingleHistoRRF (public)
+//--------------------------------------------------------------------------
+/**
+ * <p>Get a processed single histogram RRF data set.
+ *
+ * <b>return:</b>
+ * - pointer to the run data set (processed data) if data set is found
+ * - null pointer otherwise
+ *
+ * \param index msr-file run index
+ * \param tag kIndex -> data at index, kRunNo -> data of given run no
+ */
+PRunData* PRunListCollection::GetSingleHistoRRF(UInt_t index, EDataSwitch tag)
+{
+  PRunData *data = 0;
+
+  switch (tag) {
+    case kIndex:
+      if ((index < 0) || (index >= fRunSingleHistoRRFList.size())) {
+        cerr << endl << "PRunListCollection::GetSingleHistoRRF: **ERROR** index = " << index << " out of bounds";
+        cerr << endl;
+        return 0;
+      }
+
+      fRunSingleHistoRRFList[index]->CalcTheory();
+      data = fRunSingleHistoRRFList[index]->GetData();
+      break;
+    case kRunNo:
+      for (UInt_t i=0; i<fRunSingleHistoRRFList.size(); i++) {
+        if (fRunSingleHistoRRFList[i]->GetRunNo() == index) {
+          data = fRunSingleHistoRRFList[i]->GetData();
+          break;
+        }
+      }
+      break;
+    default: // error
+      break;
+  }
+
+  return data;
+}
+
 //--------------------------------------------------------------------------
 // GetAsymmetry (public)
 //--------------------------------------------------------------------------