added correct handling of alltogether unused variables

src/ToDo.txt

@@ -31,13 +31,27 @@ short term:
 
 * fix problems in FUNCTIONS of the form (sin(par1)) **DONE** 08-02-18
 
-* fix output problems of the mlog files
+* fix output problems of the mlog files **DONE** 08-02-22
-
-* implement table based theory functions (LF stuff)
 
 * if a parameter is not used at all, minuit is still varying it!! This is stupid.
   Check the minuit manual, there must be a function to "remove" these parameters,
-  i.e. forcing minuit to ignore them.
+  i.e. forcing minuit to ignore them. checked the manual: use fix() (see p.31)
+  **DONE** 08-03-10
+
+* implement table based theory functions (LF stuff)
+
+* check midas keyboard routines for the usability
+
+  system.c
+  while (ss_kbhit()) {
+      ch = ss_getchar(0);
+      if (ch == -1)
+        ch = getchar();
+
+      if (ch == '!')
+        status = RPC_SHUTDOWN;
+  }
+
 
 ---------------------
 intermediate term:

src/classes/PFitter.cpp

@@ -250,6 +250,14 @@ bool PFitter::SetParameters()
     }
   }
 
+  // check if there is an unused parameter, if so, fix it
+  for (unsigned int i=0; i<fParams.size(); i++) {
+    if (fRunInfo->ParameterInUse(i) == 0) { // parameter not used in the whole theory!!
+      fMnUserParams.Fix(i); // fix the unused parameter so that minuit will not vary it
+      cout << endl << "**WARNING** : Parameter No " << i+1 << " is not used at all, will fix it" << endl;
+    }
+  }
+
   return true;
 }
 
@@ -391,7 +399,9 @@ bool PFitter::ExecuteMinos()
 
   for (unsigned int i=0; i<fParams.size(); i++) {
     // only try to call minos if the parameter is not fixed!!
-    if (fMnUserParams.Error(i) != 0) {
+    // the 1st condition is from an user fixed variable,
+    // the 2nd condition is from an all together unused variable
+    if ((fMnUserParams.Error(i) != 0) && (fRunInfo->ParameterInUse(i) != 0)) {
       // 1-sigma MINOS errors
       std::pair<double, double> err = minos(i);
 

src/classes/PMsrHandler.cpp

@@ -82,6 +82,7 @@ PMsrHandler::~PMsrHandler()
   fCommands.clear();
   fPlots.clear();
   fStatistic.fStatLines.clear();
+  fParamInUse.clear();
 }
 
 //--------------------------------------------------------------------------
@@ -208,6 +209,9 @@ int PMsrHandler::ReadMsrFile()
   if (!HandleStatisticEntry(statistic))
     error = PMUSR_MSR_SYNTAX_ERROR;
 
+  // fill parameter in use vector
+  FillParameterInUse(theory, functions, run);
+
   // clean up
   fit_parameter.clear();
   theory.clear();
@@ -694,6 +698,26 @@ bool PMsrHandler::SetMsrParamPosError(unsigned int i, double value)
   return true;
 }
 
+//--------------------------------------------------------------------------
+// ParameterInUse (public)
+//--------------------------------------------------------------------------
+/**
+ * <p>Needed for the following purpose: if minuit is minimizing, it varies
+ * all the parameters of the parameter list (if not fixed), even if a particular
+ * parameter is <b>NOT</b> used at all. This is stupid! Hence one has to check
+ * if the parameter is used at all and if not, it has to be fixed.
+ *
+ * \param paramNo
+ */
+int PMsrHandler::ParameterInUse(unsigned int paramNo)
+{
+  // check that paramNo is within acceptable range
+  if ((paramNo < 0) || (paramNo > fParam.size()))
+    return -1;
+
+  return fParamInUse[paramNo];
+}
+
 //--------------------------------------------------------------------------
 // HandleFitParameterEntry (private)
 //--------------------------------------------------------------------------
@@ -1779,4 +1803,304 @@ bool PMsrHandler::HandleStatisticEntry(PMsrLines &lines)
   return true;
 }
 
+//--------------------------------------------------------------------------
+// FillParameterInUse (private)
+//--------------------------------------------------------------------------
+/**
+ * <p>
+ *
+ * \param theory
+ * \param funcs
+ * \param run
+ */
+void PMsrHandler::FillParameterInUse(PMsrLines &theory, PMsrLines &funcs, PMsrLines &run)
+{
+  PIntVector map;
+  PIntVector fun;
+  PMsrLines::iterator iter;
+  TObjArray  *tokens = 0;
+  TObjString *ostr = 0;
+  TString    str;
+  int        ival;
+
+  // create and initialize fParamInUse vector
+  for (unsigned int i=0; i<fParam.size(); i++)
+    fParamInUse.push_back(0);
+
+  // go through all the theory lines ------------------------------------
+  for (iter = theory.begin(); iter != theory.end(); ++iter) {
+    // remove potential comments
+    str = iter->fLine;
+    if (str.First('#') != -1)
+      str.Resize(str.First('#'));
+
+    // everything to lower case
+    str.ToLower();
+
+    // tokenize string
+    tokens = str.Tokenize(" \t");
+    if (!tokens)
+      continue;
+
+    // filter param no, map no, and fun no
+    for (int i=0; i<tokens->GetEntries(); i++) {
+      ostr = dynamic_cast<TObjString*>(tokens->At(i));
+      str = ostr->GetString();
+      if (str.IsDigit()) { // parameter number
+        ival = str.Atoi();
+        if ((ival > 0) && (ival < (int)fParam.size())) {
+          fParamInUse[ival-1]++;
+//cout << endl << ">>>> theo: param no : " << ival;
+        }
+      } else if (str.Contains("map")) { // map
+        if (FilterFunMapNumber(str, "map", ival))
+          map.push_back(ival-MSR_PARAM_MAP_OFFSET);
+      } else if (str.Contains("fun")) { // fun
+        if (FilterFunMapNumber(str, "fun", ival))
+          fun.push_back(ival-MSR_PARAM_FUN_OFFSET);
+      }
+    }
+
+    // delete tokens
+    if (tokens) {
+      delete tokens;
+      tokens = 0;
+    }
+  }
+
+  // go through all the function lines: 1st time -----------------------------
+  for (iter = funcs.begin(); iter != funcs.end(); ++iter) {
+    // remove potential comments
+    str = iter->fLine;
+    if (str.First('#') != -1)
+      str.Resize(str.First('#'));
+
+    // everything to lower case
+    str.ToLower();
+
+    tokens = str.Tokenize(" /t");
+    if (!tokens)
+      continue;
+
+    // filter fun number
+    ostr = dynamic_cast<TObjString*>(tokens->At(0));
+    str = ostr->GetString();
+    if (!FilterFunMapNumber(str, "fun", ival))
+      continue;
+    ival -= MSR_PARAM_FUN_OFFSET;
+
+    // check if fun number is used, and if yes, filter parameter numbers and maps
+    TString sstr;
+    for (unsigned int i=0; i<fun.size(); i++) {
+      if (fun[i] == ival) { // function number found
+        // filter for parX
+        sstr = iter->fLine;
+        char sval[128];
+        while (sstr.Index("par") != -1) {
+          memset(sval, 0, sizeof(sval));
+          sstr = &sstr[sstr.Index("par")+3]; // trunc sstr
+          for (int j=0; j<sstr.Sizeof(); j++) {
+            if (!isdigit(sstr[j]))
+              break;
+            sval[j] = sstr[j];
+          }
+          sscanf(sval, "%d", &ival);
+//cout << endl << ">>>> parX from func 1st, X = " << ival;
+          fParamInUse[ival-1]++;
+        }
+
+        // filter for mapX
+        sstr = iter->fLine;
+        while (sstr.Index("map") != -1) {
+          memset(sval, 0, sizeof(sval));
+          sstr = &sstr[sstr.Index("map")+3]; // trunc sstr
+          for (int j=0; j<sstr.Sizeof(); j++) {
+            if (!isdigit(sstr[j]))
+              break;
+            sval[j] = sstr[j];
+          }
+          sscanf(sval, "%d", &ival);
+//cout << endl << ">>>> mapX from func 1st, X = " << ival;
+          // check if map value already in map, otherwise add it
+          if (ival > 0) {
+            unsigned int pos;
+            for (pos=0; pos<map.size(); pos++) {
+              if (ival == map[pos])
+                break;
+            }
+            if ((unsigned int)ival == map.size()) { // new map value
+              map.push_back(ival);
+            }
+          }
+        }
+      }
+    }
+
+    // delete tokens
+    if (tokens) {
+      delete tokens;
+      tokens = 0;
+    }
+  }
+
+  // go through all the run block lines -------------------------------------
+  for (iter = run.begin(); iter != run.end(); ++iter) {
+    // remove potential comments
+    str = iter->fLine;
+    if (str.First('#') != -1)
+      str.Resize(str.First('#'));
+
+    // everything to lower case
+    str.ToLower();
+
+    // handle everything but the maps
+    if (str.Contains("alpha") || str.Contains("beta") ||
+        str.Contains("alpha2") || str.Contains("beta2") ||
+        str.Contains("norm") || str.Contains("backgr.fit") ||
+        str.Contains("rphase") || str.Contains("lifetime ")) {
+      // tokenize string
+      tokens = str.Tokenize(" \t");
+      if (!tokens)
+        continue;
+      if (tokens->GetEntries()<2)
+        continue;
+
+      ostr = dynamic_cast<TObjString*>(tokens->At(1)); // parameter number or function
+      str = ostr->GetString();
+      // check if parameter number
+      if (str.IsDigit()) {
+        ival = str.Atoi();
+        fParamInUse[ival-1]++;
+//cout << endl << ">>>> run : parameter no : " << ival;
+      }
+      // check if fun
+      if (str.Contains("fun")) {
+        if (FilterFunMapNumber(str, "fun", ival)) {
+          fun.push_back(ival-MSR_PARAM_FUN_OFFSET);
+//cout << endl << ">>>> run : fun no : " << ival-MSR_PARAM_FUN_OFFSET;
+        }
+      }
+
+      // delete tokens
+      if (tokens) {
+        delete tokens;
+        tokens = 0;
+      }
+    }
+
+    // handle the maps
+    if (str.Contains("map")) {
+      // tokenize string
+      tokens = str.Tokenize(" \t");
+      if (!tokens)
+        continue;
+
+      // get the parameter number via map
+      for (unsigned int i=0; i<map.size(); i++) {
+        if (map[i] < tokens->GetEntries()) {
+          ostr = dynamic_cast<TObjString*>(tokens->At(map[i]));
+          str = ostr->GetString();
+          if (str.IsDigit()) {
+            ival = str.Atoi();
+            fParamInUse[ival-1]++; // this is OK since map is ranging from 1 ..
+//cout << endl << ">>>> param no : " << ival << ", via map no : " << map[i];
+          }
+        }
+      }
+
+      // delete tokens
+      if (tokens) {
+        delete tokens;
+        tokens = 0;
+      }
+    }
+  }
+
+  // go through all the function lines: 2nd time -----------------------------
+  for (iter = funcs.begin(); iter != funcs.end(); ++iter) {
+    // remove potential comments
+    str = iter->fLine;
+    if (str.First('#') != -1)
+      str.Resize(str.First('#'));
+
+    // everything to lower case
+    str.ToLower();
+
+    tokens = str.Tokenize(" /t");
+    if (!tokens)
+      continue;
+
+    // filter fun number
+    ostr = dynamic_cast<TObjString*>(tokens->At(0));
+    str = ostr->GetString();
+    if (!FilterFunMapNumber(str, "fun", ival))
+      continue;
+    ival -= MSR_PARAM_FUN_OFFSET;
+
+    // check if fun number is used, and if yes, filter parameter numbers and maps
+    TString sstr;
+    for (unsigned int i=0; i<fun.size(); i++) {
+      if (fun[i] == ival) { // function number found
+        // filter for parX
+        sstr = iter->fLine;
+        char sval[128];
+        while (sstr.Index("par") != -1) {
+          memset(sval, 0, sizeof(sval));
+          sstr = &sstr[sstr.Index("par")+3]; // trunc sstr
+          for (int j=0; j<sstr.Sizeof(); j++) {
+            if (!isdigit(sstr[j]))
+              break;
+            sval[j] = sstr[j];
+          }
+          sscanf(sval, "%d", &ival);
+//cout << endl << ">>>> parX from func 2nd, X = " << ival;
+          fParamInUse[ival-1]++;
+        }
+
+        // filter for mapX
+        sstr = iter->fLine;
+        while (sstr.Index("map") != -1) {
+          memset(sval, 0, sizeof(sval));
+          sstr = &sstr[sstr.Index("map")+3]; // trunc sstr
+          for (int j=0; j<sstr.Sizeof(); j++) {
+            if (!isdigit(sstr[j]))
+              break;
+            sval[j] = sstr[j];
+          }
+          sscanf(sval, "%d", &ival);
+//cout << endl << ">>>> mapX from func 2nd, X = " << ival;
+          // check if map value already in map, otherwise add it
+          if (ival > 0) {
+            unsigned int pos;
+            for (pos=0; pos<map.size(); pos++) {
+              if (ival == map[pos])
+                break;
+            }
+            if ((unsigned int)ival == map.size()) { // new map value
+              map.push_back(ival);
+            }
+          }
+        }
+      }
+    }
+
+    // delete tokens
+    if (tokens) {
+      delete tokens;
+      tokens = 0;
+    }
+  }
+
+/*
+cout << endl << ">> fParamInUse: ";
+for (unsigned int i=0; i<fParamInUse.size(); i++)
+  cout << fParamInUse[i] << " ";
+cout << endl;
+*/
+
+  // clean up
+  map.clear();
+  fun.clear();
+}
+
 // end ---------------------------------------------------------------------

src/classes/PTheory.cpp

@@ -1008,20 +1008,20 @@ double PTheory::SkewedGauss(register double t, const vector<double>& paramValues
 
 
   // val[2] = sigma-, val[3] = sigma+
-  double z1 = val[3]*t/SQRT_TWO; // sigma+
+  double zp = val[3]*t/SQRT_TWO; // sigma+
-  double z2 = val[2]*t/SQRT_TWO; // sigma-
+  double zm = val[2]*t/SQRT_TWO; // sigma-
-  double g1 = TMath::Exp(-0.5*TMath::Power(t*val[3], 2.0)); // gauss sigma+
+  double gp = TMath::Exp(-0.5*TMath::Power(t*val[3], 2.0)); // gauss sigma+
-  double g2 = TMath::Exp(-0.5*TMath::Power(t*val[2], 2.0)); // gauss sigma-
+  double gm = TMath::Exp(-0.5*TMath::Power(t*val[2], 2.0)); // gauss sigma-
+  double wp = val[3]/(val[2]+val[3]); // sigma+ / (sigma+ + sigma-)
 
-  if ((z1 >= 25.0) || (z2 >= 25.0)) // needed to prevent crash of 1F1
+  if ((zp >= 25.0) || (zm >= 25.0)) // needed to prevent crash of 1F1
     skg = 2.0e300;
   else
-    skg = 0.5*TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) * ( g1 + g2 ) -
+    skg = (1.0-wp) * (gm * (TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) +
-          0.5*TMath::Sin(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) *
+                            TMath::Sin(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) * 2.0*zm/SQRT_PI*gsl_sf_hyperg_1F1(0.5,1.5,zm*zm)))
-          (
+          +
-          g1*2.0*z1/SQRT_PI*gsl_sf_hyperg_1F1(0.5,1.5,z1*z1) -
+          wp *       (gp * (TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) -
-          g2*2.0*z2/SQRT_PI*gsl_sf_hyperg_1F1(0.5,1.5,z2*z2)
+                            TMath::Sin(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) * 2.0*zp/SQRT_PI*gsl_sf_hyperg_1F1(0.5,1.5,zp*zp)));
-          );
 
   return skg;
 }

src/include/PMsrHandler.h

@@ -84,6 +84,8 @@ class PMsrHandler
     virtual double EvalFunc(unsigned int i, vector<int> map, vector<double> param)
                        { return fFuncHandler->Eval(i,map,param); }
 
+    virtual int ParameterInUse(unsigned int paramNo);
+
   private:
     TString                fFileName;
     TString                fTitle;     ///< holds the title string of the msr-file
@@ -99,6 +101,8 @@ class PMsrHandler
 
     PFunctionHandler *fFuncHandler; ///< needed to parse functions
 
+    PIntVector fParamInUse; ///< array holding the information if a particular parameter is used at all, i.e. if the theory is using it (perhaps via maps or functions)
+
     virtual bool HandleFitParameterEntry(PMsrLines &line);
     virtual bool HandleTheoryEntry(PMsrLines &line);
     virtual bool HandleFunctionsEntry(PMsrLines &line);
@@ -107,6 +111,8 @@ class PMsrHandler
     virtual bool HandlePlotEntry(PMsrLines &line);
     virtual bool HandleStatisticEntry(PMsrLines &line);
 
+    virtual void FillParameterInUse(PMsrLines &theory, PMsrLines &funcs, PMsrLines &run);
+
     virtual void InitRunParameterStructure(PMsrRunStructure &param);
     virtual bool FilterFunMapNumber(TString str, const char *filter, int &no);
 };

src/tests/skewedGaussianTest/fakeData.cpp

@@ -51,12 +51,12 @@ typedef vector<int>    PIntVector;
 
 void fakeDataSyntax()
 {
-  cout << endl << "usage: fakeData <pB-File> <paramInputFile> <outputFile> [--dump <dumpFln>]";
+  cout << endl << "usage: fakeData <pB-File> <paramInputFile> <outputFile> [--nemu] [--dump <dumpFln>]";
   cout << endl << "  <pB-File>:        file holding p(B) in columns B, pB, separated by ','";
   cout << endl << "                    comment lines start with a '#'";
   cout << endl << "  <paramInputFile>: holding parameters needed to generate the histograms";
   cout << endl << "                    (see below)";
-  cout << endl << "  <rootOutputFile>: output file name of the fake data histo file";
+  cout << endl << "  <outputFile>    : output file name of the fake data histo file";
   cout << endl << "  --dump <dumpFln>: option to dump p(B) and A_i(t) in a separate root file";
   cout << endl;
   cout << endl << "<paramInputFile> structure:";