Merged muonspin/musrfit:root6 into master

CMakeLists.txt

@@ -5,7 +5,7 @@ if (CMAKE_VERSION GREATER_EQUAL 3.12)
   cmake_policy(SET CMP0075 NEW)
 endif (CMAKE_VERSION GREATER_EQUAL 3.12)
 
-project(musrfit VERSION 1.7.4 LANGUAGES C CXX)
+project(musrfit VERSION 1.7.5 LANGUAGES C CXX)
 
 #--- musrfit specific options -------------------------------------------------
 option(nexus "build optional NeXus support. Needed for ISIS" OFF)

src/classes/PFitter.cpp

@@ -313,6 +313,9 @@ PFitter::PFitter(PMsrHandler *runInfo, PRunListCollection *runListCollection, Bo
     return;
   }
 
+  // create phase bool array
+  GetPhaseParams();
+
   // create fit function object
   fFitterFcn = new PFitterFcn(runListCollection, fUseChi2);
   if (!fFitterFcn) {
@@ -346,6 +349,209 @@ PFitter::~PFitter()
   }
 }
 
+//--------------------------------------------------------------------------
+// GetPhaseParams (private)
+//--------------------------------------------------------------------------
+/**
+ * <p>Checks which parameters are phases. This information is needed to
+ * restrict the phases to the intervall -360 to +360 degrees.
+ */
+void PFitter::GetPhaseParams()
+{
+  fPhase.resize(fRunInfo->GetNoOfParams());
+  for (unsigned int i=0; i<fPhase.size(); i++)
+    fPhase[i] = false;
+
+  // analyze theory block for parameters. Phases are present in the following
+  // default functions:
+  // user functions cannot be checked!
+  PMsrLines *theo = fRunInfo->GetMsrTheory();
+  TObjArray *tok = nullptr;
+  TObjString *ostr = nullptr;
+  TString str;
+  int pos = -1;
+  for (unsigned int i=0; i<theo->size(); i++) {
+    pos = -1;
+    TString line = theo->at(i).fLine;
+    if (line.Contains("TFieldCos") || line.Contains("tf ") ||
+        line.Contains("bessel") || line.Contains("b ") ||
+        line.Contains("skewedGss") || line.Contains("skg ") ||
+        line.Contains("staticNKTF") || line.Contains("snktf ") ||
+        line.Contains("dynamicNKTF") || line.Contains("dnktf ")) { // phase is 1st param
+      pos = 1;
+    }
+    if (line.Contains("internFld") || line.Contains("if ") ||
+        line.Contains("internBsl") || line.Contains("ib ")) { // phase is 2nd param
+      pos = 2;
+    }
+    if (line.Contains("muMinusExpTF") || line.Contains("mmsetf ")) { // phase is 5th param
+      pos = 5;
+    }
+
+    if (pos == -1)
+      continue;
+
+    // extract phase token
+    tok = line.Tokenize(" \t");
+    if (tok == nullptr) {
+      std::cerr << "PFitter::GetPhaseParams(): **ERROR** couldn't tokenize theory line string." << std::endl;
+      return;
+    }
+    if (tok->GetEntries() > pos) {
+      ostr = dynamic_cast<TObjString*>(tok->At(pos));
+      str = ostr->GetString();
+    }
+    // clean up
+    delete tok;
+    tok = nullptr;
+
+    // decode phase token. It can be funX, mapX, or a number
+    if (str.Contains("fun")) { // function
+      PIntVector parVec = GetParFromFun(str);
+      for (int i=0; i<parVec.size(); i++) {
+        if (parVec[i] <= fRunInfo->GetNoOfParams())
+          fPhase[parVec[i]-1] = true;
+      }
+    } else if (str.Contains("map")) { // map
+      PIntVector parVec = GetParFromMap(str);
+      for (int i=0; i<parVec.size(); i++) {
+        if (parVec[i] <= fRunInfo->GetNoOfParams())
+          fPhase[parVec[i]-1] = true;
+      }
+    } else { // must be a number
+      int idx = str.Atoi();
+      if (idx == 0) { // something went wrong, str is not an integer
+        std::cerr << "PFitter::GetPhaseParams(): **ERROR** str=" << str.View() << " is not an integer!" << std::endl;
+        return;
+      }
+      idx -= 1; // param start at 1, vector at 0
+      if (idx >= fRunInfo->GetNoOfParams()) { // idx is out-of-range
+        std::cerr << "PFitter::GetPhaseParams(): **ERROR** idx=" << idx << " is > #param = " << fRunInfo->GetNoOfParams() << "!" << std::endl;
+        return;
+      }
+      fPhase[idx] = true;
+    }
+  }
+}
+
+//--------------------------------------------------------------------------
+// GetParFromFun (private)
+//--------------------------------------------------------------------------
+/**
+ * <p>Extract from string funX the function number. Base on the function number
+ * the paramter numbers will be collected.
+ *
+ * @param funStr string of the form funX, where X is the function number
+ * @return a vector of all the parameter numbers related to funX
+ */
+PIntVector PFitter::GetParFromFun(const TString funStr)
+{
+  PIntVector parVec;
+
+  PMsrLines *funList = fRunInfo->GetMsrFunctions();
+  TObjArray *tok = nullptr;
+  TObjString *ostr = nullptr;
+  TString str;
+
+  for (int i=0; i<funList->size(); i++) {
+    if (funList->at(i).fLine.Contains(funStr)) {
+      // tokenize function string
+      tok = funList->at(i).fLine.Tokenize(" =+-*/");
+      if (tok == nullptr) {
+        std::cerr << "PFitter::GetParFromFun(): **ERROR** couldn't tokenize function string." << std::endl;
+        return parVec;
+      }
+
+      for (int j=1; j<tok->GetEntries(); j++) {
+        ostr = dynamic_cast<TObjString*>(tok->At(j));
+        str = ostr->GetString();
+        // parse tok for parX
+        if (str.Contains("par")) {
+          // find start idx of par in token
+          Ssiz_t idx = str.Index("par");
+          idx += 3;
+          TString parStr("");
+          do {
+            parStr += str[idx];
+          } while (isdigit(str[idx++]));
+          parVec.push_back(parStr.Atoi());
+        }
+        // parse tok for mapX
+        if (str.Contains("map")) {
+          // find start idx of par in token
+          Ssiz_t idx = str.Index("map");
+          idx += 3;
+          TString mapStr("map");
+          do {
+            mapStr += str[idx];
+          } while (isdigit(str[idx++]));
+          PIntVector mapParVec = GetParFromMap(mapStr);
+          for (int k=0; k<mapParVec.size(); k++) {
+            parVec.push_back(mapParVec[k]);
+          }
+        }
+      }
+
+      // clean up
+      delete tok;
+      tok = nullptr;
+    }
+  }
+
+  return parVec;
+}
+
+//--------------------------------------------------------------------------
+// GetParFromMap (private)
+//--------------------------------------------------------------------------
+/**
+ * <p>Extract from string mapX the map number. Based on the map number the
+ * parameter numbers will be collected.
+ *
+ * @param mapStr string of the form mapX, where X is the map number
+ * @return a vector of all the parameter numbers related to mapX
+ */
+PIntVector PFitter::GetParFromMap(const TString mapStr)
+{
+  PIntVector parVec;
+
+  TString str = mapStr;
+  str.Remove(0,3); // remove map from string
+
+  int idx=str.Atoi();
+  if (idx == 0) {
+    std::cerr << "PFitter::GetParFromMap(): **ERROR** couldn't get propper index from mapX!" << std::endl;
+    return parVec;
+  }
+
+  idx -= 1; // map starts at 1, map vector at 0
+
+  // go through all the runs and collect the parameters from the map vectors
+  PMsrRunList *runList = fRunInfo->GetMsrRunList();
+  if (runList == nullptr) {
+    std::cerr << "PFitter::GetParFromMap(): **ERROR** couldn't get required run list information!" << std::endl;
+    return parVec;
+  }
+
+  PIntVector *map = nullptr;
+  for (int i=0; i<runList->size(); i++) {
+    map = runList->at(i).GetMap();
+    if (map == nullptr) {
+      std::cerr << "PFitter::GetParFromMap(): **ERROR** couldn't get required map information (idx=" << i << ")!" << std::endl;
+      parVec.clear();
+      return parVec;
+    }
+    if (idx >= map->size()) {
+      std::cerr << "PFitter::GetParFromMap(): **ERROR** requested map index (idx=" << idx << ") out-of-range (" << map->size() << ")!" << std::endl;
+      parVec.clear();
+      return parVec;
+    }
+    parVec.push_back(map->at(idx));
+  }
+
+  return parVec;
+}
+
 //--------------------------------------------------------------------------
 // DoFit
 //--------------------------------------------------------------------------
@@ -1608,7 +1814,12 @@ Bool_t PFitter::ExecuteMigrad()
 
   // fill run info
   for (UInt_t i=0; i<fParams.size(); i++) {
-    fRunInfo->SetMsrParamValue(i, min.UserState().Value(i));
+    Double_t dval = min.UserState().Value(i);
+    if (fPhase[i]) {
+      Int_t m = (Int_t)(dval/360.0);
+      dval = dval - m*360.0;
+    }
+    fRunInfo->SetMsrParamValue(i, dval);
     fRunInfo->SetMsrParamStep(i, min.UserState().Error(i));
     fRunInfo->SetMsrParamPosErrorPresent(i, false);
   }
@@ -1684,7 +1895,12 @@ Bool_t PFitter::ExecuteMinimize()
 
   // fill run info
   for (UInt_t i=0; i<fParams.size(); i++) {
-    fRunInfo->SetMsrParamValue(i, min.UserState().Value(i));
+    Double_t dval = min.UserState().Value(i);
+    if (fPhase[i]) {
+      Int_t m = (Int_t)(dval/360.0);
+      dval = dval - m*360.0;
+    }
+    fRunInfo->SetMsrParamValue(i, dval);
     fRunInfo->SetMsrParamStep(i, min.UserState().Error(i));
     fRunInfo->SetMsrParamPosErrorPresent(i, false);
   }
@@ -2392,7 +2608,12 @@ Bool_t PFitter::ExecuteSimplex()
 
   // fill run info
   for (UInt_t i=0; i<fParams.size(); i++) {
-    fRunInfo->SetMsrParamValue(i, min.UserState().Value(i));
+    Double_t dval = min.UserState().Value(i);
+    if (fPhase[i]) {
+      Int_t m = (Int_t)(dval/360.0);
+      dval = dval - m*360.0;
+    }
+    fRunInfo->SetMsrParamValue(i, dval);
     fRunInfo->SetMsrParamStep(i, min.UserState().Error(i));
     fRunInfo->SetMsrParamPosErrorPresent(i, false);
   }

src/include/PFitter.h

@@ -155,6 +155,13 @@ class PFitter
     Bool_t fSectorFlag; ///< sector command present flag    
     std::vector<PSectorChisq> fSector; ///< stores all chisq/maxLH sector information
 
+    std::vector<bool> fPhase; ///< flag array in which an entry is true if the related parameter value is a phase
+
+    // phase related functions
+    void GetPhaseParams();
+    PIntVector GetParFromFun(const TString funStr);
+    PIntVector GetParFromMap(const TString mapStr);
+
     // commands
     Bool_t CheckCommands();
     Bool_t SetParameters();