diff --git a/CMakeLists.txt b/CMakeLists.txt index aef7326a..3875fe16 100644 --- a/CMakeLists.txt +++ b/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.6.0 LANGUAGES C CXX) +project(musrfit VERSION 1.6.1 LANGUAGES C CXX) #--- musrfit specific options ------------------------------------------------- option(nexus "build optional NeXus support. Needed for ISIS" OFF) diff --git a/ChangeLog b/ChangeLog index 83f6f199..074bbde2 100644 --- a/ChangeLog +++ b/ChangeLog @@ -12,6 +12,11 @@ or https://bitbucket.org/muonspin/musrfit/commits/all +Release of V1.6.1, 2020/05/25 +============================= + +For asymmetry fits: allow alpha, beta in the RUN block to be functions. + Release of V1.6.0, 2020/05/16 ============================= diff --git a/doc/html/.buildinfo b/doc/html/.buildinfo index 7ac2c58d..b66b59df 100644 --- a/doc/html/.buildinfo +++ b/doc/html/.buildinfo @@ -1,4 +1,4 @@ # Sphinx build info version 1 # This file hashes the configuration used when building these files. When it is not found, a full rebuild will be done. -config: 90c2680fff90669033f1bce5f6912272 +config: 2fd2e79340d6214c8fde0515f4f3fbe3 tags: 645f666f9bcd5a90fca523b33c5a78b7 diff --git a/doc/html/_sources/user-manual.rst.txt b/doc/html/_sources/user-manual.rst.txt index 3e8e23bc..b20c63dd 100644 --- a/doc/html/_sources/user-manual.rst.txt +++ b/doc/html/_sources/user-manual.rst.txt @@ -775,7 +775,7 @@ The FUNCTIONS Block +++++++++++++++++++ Here some auxiliary functions can be defined. These functions can currently *only* operate on the defined parameters. They can be used in the -:ref:`THEORY block ` and for one specific case in the :ref:`RUN block `. Supported is the use of basic arithmetic: +:ref:`THEORY block ` and for three specific cases in the :ref:`RUN block ` (`norm`, `alpha`, and `beta`). Supported is the use of basic arithmetic: :math:`+` Addition @@ -1064,18 +1064,24 @@ In order to describe the operations needed for fitting and plotting, quite some .. _msr-alpha-beta: **alpha, beta** (fit type 2, 3, 5) - These parameters are used to correct the asymmetry for different detector efficiencies, solid angles and initial asymmetries. They are defined as :math:`\alpha = N_{0,b}/N_{0,f}` and :math:`\beta = A_{0,b}/A_{0,f}`. If the parameters are not specified in the :ref:`RUN block `, for each one the value of 1 is assumed (for fittype 5, alpha is estimated from the ratio of :math:`\sum_i \left( N_{\mathrm{bp}}(i)+N_{\mathrm{bm}}(i) \right)` and :math:`\sum_i \left( N_{\mathrm{fp}}(i)+N_{\mathrm{fm}}(i) \right)`). Example for alpha with fit parameter number 1: + These parameters are used to correct the asymmetry for different detector efficiencies, solid angles and initial asymmetries. They are defined as :math:`\alpha = N_{0,b}/N_{0,f}` and :math:`\beta = A_{0,b}/A_{0,f}`. If the parameters are not specified in the :ref:`RUN block `, for each one the value of 1 is assumed (for fittype 5, alpha is estimated from the ratio of :math:`\sum_i \left( N_{\mathrm{bp}}(i)+N_{\mathrm{bm}}(i) \right)` and :math:`\sum_i \left( N_{\mathrm{fp}}(i)+N_{\mathrm{fm}}(i) \right)`). Both, `alpha` as well as `beta` can be expressed through a function. Example for alpha with fit parameter number 1: :: alpha 1 + Example for an ``alpha`` defined via function number 1: + + :: + + alpha fun1 + .. index:: norm .. _msr-norm: **norm** (fit type 0) Number of the fit parameter that represents the normalization constant :math:`N_0` of the histogram; the value of this parameter is given either per nanosecond or per bin (see :ref:`below `). - It is possible to substitute the parameter number by a function here (and only here in a RUN block), for instance to relate :math:`N_0`\'s of different histograms through an :math:`\alpha` + It is possible to substitute the parameter number by a function, for instance to relate :math:`N_0`\'s of different histograms through an :math:`\alpha` parameter. Example for a ``norm`` defined by fit parameter number 12: :: diff --git a/doc/html/_static/documentation_options.js b/doc/html/_static/documentation_options.js index c5d9c26f..63761070 100644 --- a/doc/html/_static/documentation_options.js +++ b/doc/html/_static/documentation_options.js @@ -1,6 +1,6 @@ var DOCUMENTATION_OPTIONS = { URL_ROOT: document.getElementById("documentation_options").getAttribute('data-url_root'), - VERSION: '1.6.0', + VERSION: '1.6.1', LANGUAGE: 'None', COLLAPSE_INDEX: false, FILE_SUFFIX: '.html', diff --git a/doc/html/acknowledgement.html b/doc/html/acknowledgement.html index cfb036e7..ffebb3a1 100644 --- a/doc/html/acknowledgement.html +++ b/doc/html/acknowledgement.html @@ -4,7 +4,7 @@ - Acknowledgements — musrfit 1.6.0 documentation + Acknowledgements — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -100,12 +100,12 @@ extremely competent way to deal with his projects as well as to deal with the ch
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/any2many.html b/doc/html/any2many.html index bc76f097..ac86b1bd 100644 --- a/doc/html/any2many.html +++ b/doc/html/any2many.html @@ -4,7 +4,7 @@ - any2many - a Universal μSR-file-format converter — musrfit 1.6.0 documentation + any2many - a Universal μSR-file-format converter — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -90,12 +90,12 @@ For a detailed description see previous | -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/bugtracking.html b/doc/html/bugtracking.html index d0421f76..e13f9a74 100644 --- a/doc/html/bugtracking.html +++ b/doc/html/bugtracking.html @@ -4,7 +4,7 @@ - Bugtracking — musrfit 1.6.0 documentation + Bugtracking — musrfit 1.6.1 documentation @@ -26,7 +26,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -81,12 +81,12 @@ or send an e-mail to A. Suter at PSI.

  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/cite.html b/doc/html/cite.html index 4c1e816f..d9161bea 100644 --- a/doc/html/cite.html +++ b/doc/html/cite.html @@ -4,7 +4,7 @@ - How to Cite musrfit? — musrfit 1.6.0 documentation + How to Cite musrfit? — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -102,12 +102,12 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/genindex.html b/doc/html/genindex.html index e471789c..327beab5 100644 --- a/doc/html/genindex.html +++ b/doc/html/genindex.html @@ -5,7 +5,7 @@ - Index — musrfit 1.6.0 documentation + Index — musrfit 1.6.1 documentation @@ -23,7 +23,7 @@
  • index
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -611,12 +611,12 @@
  • index
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/index.html b/doc/html/index.html index 03856aa8..d2e26fb2 100644 --- a/doc/html/index.html +++ b/doc/html/index.html @@ -4,7 +4,7 @@ - Welcome to the musrfit documentation! — musrfit 1.6.0 documentation + Welcome to the musrfit documentation! — musrfit 1.6.1 documentation @@ -26,7 +26,7 @@
  • next |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -168,12 +168,12 @@
  • next |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/msr2data.html b/doc/html/msr2data.html index af84f1ff..0e7621e4 100644 --- a/doc/html/msr2data.html +++ b/doc/html/msr2data.html @@ -4,7 +4,7 @@ - msr2data - A Program for Automatically Processing Multiple musrfit msr Files — musrfit 1.6.0 documentation + msr2data - A Program for Automatically Processing Multiple musrfit msr Files — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -418,12 +418,12 @@ fit serves as template for the second and so on. The template field stays empty
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/mupp.html b/doc/html/mupp.html index 40de7bee..efdeaff2 100644 --- a/doc/html/mupp.html +++ b/doc/html/mupp.html @@ -4,7 +4,7 @@ - mupp - μSR Parameter Plotter — musrfit 1.6.0 documentation + mupp - μSR Parameter Plotter — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -298,12 +298,12 @@ SCRIPT COMMANDS:
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/musr-root.html b/doc/html/musr-root.html index 85684427..876998cb 100644 --- a/doc/html/musr-root.html +++ b/doc/html/musr-root.html @@ -4,7 +4,7 @@ - MusrRoot - an Extensible Open File Format for μSR — musrfit 1.6.0 documentation + MusrRoot - an Extensible Open File Format for μSR — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -910,12 +910,12 @@ the entry has been added. The last token, previous | -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/musredit.html b/doc/html/musredit.html index debed612..9352626e 100644 --- a/doc/html/musredit.html +++ b/doc/html/musredit.html @@ -4,7 +4,7 @@ - musredit: the GUI Based Interface to musrfit — musrfit 1.6.0 documentation + musredit: the GUI Based Interface to musrfit — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -531,12 +531,12 @@ the corresponding fit parameter value, except the phases where the step will be
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/search.html b/doc/html/search.html index d0cd2faa..99c1e1a2 100644 --- a/doc/html/search.html +++ b/doc/html/search.html @@ -4,7 +4,7 @@ - Search — musrfit 1.6.0 documentation + Search — musrfit 1.6.1 documentation @@ -27,7 +27,7 @@
  • index
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -75,12 +75,12 @@
  • index
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/setup-dks.html b/doc/html/setup-dks.html index c5b25e44..2532f719 100644 --- a/doc/html/setup-dks.html +++ b/doc/html/setup-dks.html @@ -4,7 +4,7 @@ - Setting up musrfit / DKS: High Speed Fitting with GPU’s — musrfit 1.6.0 documentation + Setting up musrfit / DKS: High Speed Fitting with GPU’s — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -309,12 +309,12 @@ The only thing you need previous | -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/setup-standard.html b/doc/html/setup-standard.html index 5c0f0b16..49a893e2 100644 --- a/doc/html/setup-standard.html +++ b/doc/html/setup-standard.html @@ -4,7 +4,7 @@ - Setting up musrfit on Different Platforms — musrfit 1.6.0 documentation + Setting up musrfit on Different Platforms — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -1293,12 +1293,12 @@ $ musrview test-histo-ROOT-NPP.msr
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/tutorial.html b/doc/html/tutorial.html index 8414ee1e..889b133a 100644 --- a/doc/html/tutorial.html +++ b/doc/html/tutorial.html @@ -4,7 +4,7 @@ - Tutorial for musrfit — musrfit 1.6.0 documentation + Tutorial for musrfit — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -431,12 +431,12 @@ For a complete description please refer to the manuals of previous | -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/user-libs.html b/doc/html/user-libs.html index 43df5c10..87003663 100644 --- a/doc/html/user-libs.html +++ b/doc/html/user-libs.html @@ -4,7 +4,7 @@ - Documentation of user libs (user functions) — musrfit 1.6.0 documentation + Documentation of user libs (user functions) — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -647,12 +647,12 @@ K(m)&=\int_0^{\pi/2}\frac{\mathrm d\varphi}{\sqrt{1-m^2\sin^2{\varphi}}},\en
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/doc/html/user-manual.html b/doc/html/user-manual.html index 989435e4..25bbb801 100644 --- a/doc/html/user-manual.html +++ b/doc/html/user-manual.html @@ -4,7 +4,7 @@ - User manual — musrfit 1.6.0 documentation + User manual — musrfit 1.6.1 documentation @@ -30,7 +30,7 @@
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • @@ -897,7 +897,7 @@ instead of specifying the parameters directly.

    The FUNCTIONS Block

    Here some auxiliary functions can be defined. These functions can currently only operate on the defined parameters. They can be used in the -THEORY block and for one specific case in the RUN block. Supported is the use of basic arithmetic:

    +THEORY block and for three specific cases in the RUN block (norm, alpha, and beta). Supported is the use of basic arithmetic:

    \(+\)

    Addition

    @@ -1169,15 +1169,19 @@ a short digression is needed.

    -
    alpha, beta (fit type 2, 3, 5)

    These parameters are used to correct the asymmetry for different detector efficiencies, solid angles and initial asymmetries. They are defined as \(\alpha = N_{0,b}/N_{0,f}\) and \(\beta = A_{0,b}/A_{0,f}\). If the parameters are not specified in the RUN block, for each one the value of 1 is assumed (for fittype 5, alpha is estimated from the ratio of \(\sum_i \left( N_{\mathrm{bp}}(i)+N_{\mathrm{bm}}(i) \right)\) and \(\sum_i \left( N_{\mathrm{fp}}(i)+N_{\mathrm{fm}}(i) \right)\)). Example for alpha with fit parameter number 1:

    +
    alpha, beta (fit type 2, 3, 5)

    These parameters are used to correct the asymmetry for different detector efficiencies, solid angles and initial asymmetries. They are defined as \(\alpha = N_{0,b}/N_{0,f}\) and \(\beta = A_{0,b}/A_{0,f}\). If the parameters are not specified in the RUN block, for each one the value of 1 is assumed (for fittype 5, alpha is estimated from the ratio of \(\sum_i \left( N_{\mathrm{bp}}(i)+N_{\mathrm{bm}}(i) \right)\) and \(\sum_i \left( N_{\mathrm{fp}}(i)+N_{\mathrm{fm}}(i) \right)\)). Both, alpha as well as beta can be expressed through a function. Example for alpha with fit parameter number 1:

    alpha 1
    +

    Example for an alpha defined via function number 1:

    +
    alpha fun1
+
    +
    norm (fit type 0)

    Number of the fit parameter that represents the normalization constant \(N_0\) of the histogram; the value of this parameter is given either per nanosecond or per bin (see below). -It is possible to substitute the parameter number by a function here (and only here in a RUN block), for instance to relate \(N_0\)’s of different histograms through an \(\alpha\) +It is possible to substitute the parameter number by a function, for instance to relate \(N_0\)’s of different histograms through an \(\alpha\) parameter. Example for a norm defined by fit parameter number 12:

    norm 12
    @@ -2429,12 +2433,12 @@ In case this cannot be ensured, the parallelization can be disabled by –di
  • previous |
    • -
  • musrfit 1.6.0 documentation »
    • +
  • musrfit 1.6.1 documentation »
    • diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 5bedae20..fd4f70fb 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -819,15 +819,24 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) } else if (sstr.BeginsWith("alpha ")) { fout.width(16); fout << std::left << "alpha"; - fout << fRuns[runNo].GetAlphaParamNo() << std::endl; + // check of alpha is given as a function + if (fRuns[runNo].GetAlphaParamNo() >= MSR_PARAM_FUN_OFFSET) + fout << "fun" << fRuns[runNo].GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + else + fout << fRuns[runNo].GetAlphaParamNo() << std::endl; + fout << std::endl; } else if (sstr.BeginsWith("beta ")) { fout.width(16); fout << std::left << "beta"; - fout << fRuns[runNo].GetBetaParamNo() << std::endl; + if (fRuns[runNo].GetBetaParamNo() >= MSR_PARAM_FUN_OFFSET) + fout << "fun" << fRuns[runNo].GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + else + fout << fRuns[runNo].GetBetaParamNo() << std::endl; + fout << std::endl; } else if (sstr.BeginsWith("norm")) { fout.width(16); fout << std::left << "norm"; - // check if norm is give as a function + // check if norm is given as a function if (fRuns[runNo].GetNormParamNo() >= MSR_PARAM_FUN_OFFSET) fout << "fun" << fRuns[runNo].GetNormParamNo()-MSR_PARAM_FUN_OFFSET; else @@ -1939,14 +1948,24 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, std::map= MSR_PARAM_FUN_OFFSET) + fout << "fun" << fRuns[i].GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + else + fout << fRuns[i].GetAlphaParamNo(); + fout << std::endl; } // beta if (fRuns[i].GetBetaParamNo() != -1) { fout.width(16); fout << std::left << "beta"; - fout << fRuns[i].GetBetaParamNo() << std::endl; + // check if beta is give as a function + if (fRuns[i].GetBetaParamNo() >= MSR_PARAM_FUN_OFFSET) + fout << "fun" << fRuns[i].GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + else + fout << fRuns[i].GetBetaParamNo(); + fout << std::endl; } // norm @@ -3379,6 +3398,12 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines) param.SetAlphaParamNo(ival); else error = true; + } else if (str.Contains("fun")) { + Int_t no; + if (FilterNumber(str, "fun", MSR_PARAM_FUN_OFFSET, no)) + param.SetAlphaParamNo(no); + else + error = true; } else { error = true; } @@ -3399,6 +3424,12 @@ Bool_t PMsrHandler::HandleRunEntry(PMsrLines &lines) ival = str.Atoi(); if (ival > 0) param.SetBetaParamNo(ival); + else + error = true; + } else if (str.Contains("fun")) { + Int_t no; + if (FilterNumber(str, "fun", MSR_PARAM_FUN_OFFSET, no)) + param.SetBetaParamNo(no); else error = true; } else { @@ -4582,9 +4613,12 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) // handle a single PLOT block while ((iter1 != iter2) && !error) { + TString line = iter1->fLine; + if (line.First('#') != -1) // remove trailing comment before proceed + line.Resize(line.First('#')); - if (iter1->fLine.Contains("PLOT")) { // handle plot header - tokens = iter1->fLine.Tokenize(" \t"); + if (line.Contains("PLOT")) { // handle plot header + tokens = line.Tokenize(" \t"); if (!tokens) { std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo; std::cerr << std::endl << std::endl; @@ -4605,9 +4639,9 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) delete tokens; tokens = nullptr; } - } else if (iter1->fLine.Contains("lifetimecorrection", TString::kIgnoreCase)) { + } else if (line.Contains("lifetimecorrection", TString::kIgnoreCase)) { param.fLifeTimeCorrection = true; - } else if (iter1->fLine.Contains("runs", TString::kIgnoreCase)) { // handle plot runs + } else if (line.Contains("runs", TString::kIgnoreCase)) { // handle plot runs TComplex run; PStringNumberList *rl; std::string errorMsg; @@ -4623,7 +4657,7 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) case MSR_PLOT_ASYM_RRF: case MSR_PLOT_NON_MUSR: case MSR_PLOT_MU_MINUS: - rl = new PStringNumberList(iter1->fLine.Data()); + rl = new PStringNumberList(line.Data()); if (!rl->Parse(errorMsg, true)) { std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo; std::cerr << std::endl << ">> Error Message: " << errorMsg; @@ -4643,14 +4677,14 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) error = true; break; } - } else if (iter1->fLine.Contains("range ", TString::kIgnoreCase)) { // handle plot range + } else if (line.Contains("range ", TString::kIgnoreCase)) { // handle plot range // remove previous entries param.fTmin.clear(); param.fTmax.clear(); param.fYmin.clear(); param.fYmax.clear(); - tokens = iter1->fLine.Tokenize(" \t"); + tokens = line.Tokenize(" \t"); if (!tokens) { std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo; std::cerr << std::endl << std::endl; @@ -4700,14 +4734,14 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) delete tokens; tokens = nullptr; } - } else if (iter1->fLine.Contains("sub_ranges", TString::kIgnoreCase)) { + } else if (line.Contains("sub_ranges", TString::kIgnoreCase)) { // remove previous entries param.fTmin.clear(); param.fTmax.clear(); param.fYmin.clear(); param.fYmax.clear(); - tokens = iter1->fLine.Tokenize(" \t"); + tokens = line.Tokenize(" \t"); if (!tokens) { std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo; std::cerr << std::endl << std::endl; @@ -4762,11 +4796,11 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) delete tokens; tokens = nullptr; } - } else if (iter1->fLine.Contains("use_fit_ranges", TString::kIgnoreCase)) { + } else if (line.Contains("use_fit_ranges", TString::kIgnoreCase)) { param.fUseFitRanges = true; // check if y-ranges are given - tokens = iter1->fLine.Tokenize(" \t"); + tokens = line.Tokenize(" \t"); if (!tokens) { std::cerr << std::endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize PLOT in line " << iter1->fLineNo; std::cerr << std::endl << std::endl; @@ -5207,12 +5241,20 @@ UInt_t PMsrHandler::GetNoOfFitParameters(UInt_t idx) paramVector.push_back(fRuns[idx].GetBkgFitParamNo()); // get alpha parameter if present (asymmetry fit) - if (fRuns[idx].GetAlphaParamNo() != -1) - paramVector.push_back(fRuns[idx].GetAlphaParamNo()); + if (fRuns[idx].GetAlphaParamNo() != -1) { + if (fRuns[idx].GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) // parameter + paramVector.push_back(fRuns[idx].GetAlphaParamNo()); + else // function + funVector.push_back(fRuns[idx].GetAlphaParamNo() - MSR_PARAM_FUN_OFFSET); + } // get beta parameter if present (asymmetry fit) - if (fRuns[idx].GetBetaParamNo() != -1) - paramVector.push_back(fRuns[idx].GetBetaParamNo()); + if (fRuns[idx].GetBetaParamNo() != -1) { + if (fRuns[idx].GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) // parameter + paramVector.push_back(fRuns[idx].GetBetaParamNo()); + else // function + funVector.push_back(fRuns[idx].GetBetaParamNo() - MSR_PARAM_FUN_OFFSET); + } // go through the theory block and collect parameters // possible entries: number -> parameter, fun -> function, map -> maps diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index e1586867..d4fa8e91 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -111,7 +111,16 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UIn return; } // check if alpha parameter is within proper bounds - if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > static_cast(param->size()))) { + if (fRunInfo->GetAlphaParamNo() >= MSR_PARAM_FUN_OFFSET) { // alpha seems to be a function + if ((fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET < 0) || + (fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET > msrInfo->GetNoOfFuncs())) { + std::cerr << std::endl << ">> PRunAsymmetry::PRunAsymmetry(): **ERROR** alpha parameter is a function with no = " << fRunInfo->GetAlphaParamNo(); + std::cerr << std::endl << ">> This is out of bound, since there are only " << msrInfo->GetNoOfFuncs() << " functions."; + std::cerr << std::endl; + fValid = false; + return; + } + } else if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > static_cast(param->size()))) { std::cerr << std::endl << ">> PRunAsymmetry::PRunAsymmetry(): **ERROR** alpha parameter no = " << fRunInfo->GetAlphaParamNo(); std::cerr << std::endl << ">> This is out of bound, since there are only " << param->size() << " parameters."; std::cerr << std::endl; @@ -120,10 +129,11 @@ PRunAsymmetry::PRunAsymmetry(PMsrHandler *msrInfo, PRunDataHandler *rawData, UIn } // check if alpha is fixed Bool_t alphaFixedToOne = false; - if (((*param)[fRunInfo->GetAlphaParamNo()-1].fStep == 0.0) && - ((*param)[fRunInfo->GetAlphaParamNo()-1].fValue == 1.0)) - alphaFixedToOne = true; - + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is parameter + if (((*param)[fRunInfo->GetAlphaParamNo()-1].fStep == 0.0) && + ((*param)[fRunInfo->GetAlphaParamNo()-1].fValue == 1.0)) + alphaFixedToOne = true; + } // check if beta is given Bool_t betaFixedToOne = false; if (fRunInfo->GetBetaParamNo() == -1) { // no beta given hence assuming beta == 1 @@ -196,6 +206,58 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector& par) Double_t time(1.0); Int_t i; + // determine alpha/beta + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + a = 1.0; + b = 1.0; + break; + case 2: // alpha != 1, beta == 1 + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + b = 1.0; + break; + case 3: // alpha == 1, beta != 1 + a = 1.0; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + break; + case 4: // alpha != 1, beta != 1 + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + break; + default: + a = 1.0; + b = 1.0; + break; + } + // Calculate the theory function once to ensure one function evaluation for the current set of parameters. // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once // for a given set of parameters---which should be done outside of the parallelized loop. @@ -206,34 +268,12 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector& par) Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); if (chunk < 10) chunk = 10; - #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq) + #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,f) schedule(dynamic,chunk) reduction(+:chisq) #endif for (i=fStartTimeBin; i(i)*fData.GetDataTimeStep(); - switch (fAlphaBetaTag) { - case 1: // alpha == 1, beta == 1 - asymFcnValue = fTheory->Func(time, par, fFuncValues); - break; - case 2: // alpha != 1, beta == 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); - break; - case 3: // alpha == 1, beta != 1 - b = par[fRunInfo->GetBetaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); - break; - case 4: // alpha != 1, beta != 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - b = par[fRunInfo->GetBetaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); - break; - default: - asymFcnValue = 0.0; - break; - } + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); diff = fData.GetValue()->at(i) - asymFcnValue; chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); } @@ -429,18 +469,46 @@ void PRunAsymmetry::CalcTheory() asymFcnValue = fTheory->Func(time, par, fFuncValues); break; case 2: // alpha != 1, beta == 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); break; case 3: // alpha == 1, beta != 1 - b = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); break; case 4: // alpha != 1, beta != 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - b = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); break; @@ -1138,16 +1206,44 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]) beta = 1.0; break; case 2: // alpha != 1, beta == 1 - alpha = par[fRunInfo->GetAlphaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } beta = 1.0; break; case 3: // alpha == 1, beta != 1 alpha = 1.0; - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; case 4: // alpha != 1, beta != 1 - alpha = par[fRunInfo->GetAlphaParamNo()-1]; - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; default: break; @@ -1349,16 +1445,44 @@ Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2] beta = 1.0; break; case 2: // alpha != 1, beta == 1 - alpha = par[fRunInfo->GetAlphaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } beta = 1.0; break; case 3: // alpha == 1, beta != 1 alpha = 1.0; - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; case 4: // alpha != 1, beta != 1 - alpha = par[fRunInfo->GetAlphaParamNo()-1]; - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; default: break; diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 9b95c074..66b7626c 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -206,56 +206,90 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) Double_t time(1.0),alphaest; Int_t i; + // determine alpha/beta + alphaest = fRunInfo->GetEstimatedAlpha(); + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + a = 1.0; + b = 1.0; + break; + case 2: // alpha != 1, beta == 1 + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + b = 1.0; + break; + case 3: // alpha == 1, beta != 1 + a = 1.0; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + break; + case 4: // alpha != 1, beta != 1 + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + break; + case 5: // alpha ?? , beta == 1 + a = alphaest; + b = 1.0; + break; + case 6: // alpha ??, beta != 1 + a = alphaest; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + break; + default: + a = 1.0; + b = 1.0; + break; + } + // Calculate the theory function once to ensure one function evaluation for the current set of parameters. // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once // for a given set of parameters---which should be done outside of the parallelized loop. // For all other functions it means a tiny and acceptable overhead. asymFcnValue = fTheory->Func(time, par, fFuncValues); - alphaest = fRunInfo->GetEstimatedAlpha(); #ifdef HAVE_GOMP Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); if (chunk < 10) chunk = 10; - #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq) + #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,f) schedule(dynamic,chunk) reduction(+:chisq) #endif for (i=fStartTimeBin; i(i)*fData.GetDataTimeStep(); - switch (fAlphaBetaTag) { - case 1: // alpha == 1, beta == 1 - asymFcnValue = fTheory->Func(time, par, fFuncValues); - break; - case 2: // alpha != 1, beta == 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues)/2.0; - asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0)); - break; - case 3: // alpha == 1, beta != 1 - b = par[fRunInfo->GetBetaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues)/2.0; - asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-f*(b+1.0)/(2.0+f*(b-1.0)); - break; - case 4: // alpha != 1, beta != 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - b = par[fRunInfo->GetBetaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues)/2.0; - asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); - break; - case 5: // alpha ?? , beta == 1 - a = alphaest; - f = fTheory->Func(time, par, fFuncValues)/2.0; - asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0)); - break; - case 6: // alpha ??, beta != 1 - a = alphaest; - b = par[fRunInfo->GetBetaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues)/2.0; - asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); - break; - default: - asymFcnValue = 0.0; - break; - } + f = fTheory->Func(time, par, fFuncValues)/2.0; + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); diff = fData.GetValue()->at(i) - asymFcnValue; chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); } @@ -455,18 +489,46 @@ void PRunAsymmetryBNMR::CalcTheory() asymFcnValue = fTheory->Func(time, par, fFuncValues); break; case 2: // alpha != 1, beta == 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0)); break; case 3: // alpha == 1, beta != 1 - b = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-f*(b+1.0)/(2.0+f*(b-1.0)); break; case 4: // alpha != 1, beta != 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - b = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); break; @@ -477,7 +539,14 @@ void PRunAsymmetryBNMR::CalcTheory() break; case 6: // alpha ??, beta != 1 a = alphaest; - b = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); break; @@ -1287,16 +1356,44 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 beta = 1.0; break; case 2: // alpha != 1, beta == 1 - alpha = par[fRunInfo->GetAlphaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } beta = 1.0; break; case 3: // alpha == 1, beta != 1 alpha = 1.0; - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; case 4: // alpha != 1, beta != 1 - alpha = par[fRunInfo->GetAlphaParamNo()-1]; - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; case 5: // alpha ?? , beta == 1 // use estimated value @@ -1304,7 +1401,14 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 break; case 6: // alpha ??, beta != 1 // use estimated value - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; default: break; diff --git a/src/classes/PRunAsymmetryRRF.cpp b/src/classes/PRunAsymmetryRRF.cpp index 8294c28e..4c15952b 100644 --- a/src/classes/PRunAsymmetryRRF.cpp +++ b/src/classes/PRunAsymmetryRRF.cpp @@ -188,6 +188,58 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector& par) Double_t time(1.0); Int_t i; + // determine alpha/beta + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + a = 1.0; + b = 1.0; + break; + case 2: // alpha != 1, beta == 1 + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + b = 1.0; + break; + case 3: // alpha == 1, beta != 1 + a = 1.0; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + break; + case 4: // alpha != 1, beta != 1 + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + break; + default: + a = 1.0; + b = 1.0; + break; + } + // Calculate the theory function once to ensure one function evaluation for the current set of parameters. // This is needed for the LF and user functions where some non-thread-save calculations only need to be calculated once // for a given set of parameters---which should be done outside of the parallelized loop. @@ -198,34 +250,12 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector& par) Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); if (chunk < 10) chunk = 10; - #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq) + #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,f) schedule(dynamic,chunk) reduction(+:chisq) #endif for (i=fStartTimeBin; i(i)*fData.GetDataTimeStep(); - switch (fAlphaBetaTag) { - case 1: // alpha == 1, beta == 1 - asymFcnValue = fTheory->Func(time, par, fFuncValues); - break; - case 2: // alpha != 1, beta == 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); - break; - case 3: // alpha == 1, beta != 1 - b = par[fRunInfo->GetBetaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); - break; - case 4: // alpha != 1, beta != 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - b = par[fRunInfo->GetBetaParamNo()-1]; - f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); - break; - default: - asymFcnValue = 0.0; - break; - } + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); diff = fData.GetValue()->at(i) - asymFcnValue; chisq += diff*diff / (fData.GetError()->at(i)*fData.GetError()->at(i)); } @@ -421,18 +451,46 @@ void PRunAsymmetryRRF::CalcTheory() asymFcnValue = fTheory->Func(time, par, fFuncValues); break; case 2: // alpha != 1, beta == 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); break; case 3: // alpha == 1, beta != 1 - b = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); break; case 4: // alpha != 1, beta != 1 - a = par[fRunInfo->GetAlphaParamNo()-1]; - b = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + a = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + a = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + b = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + b = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } f = fTheory->Func(time, par, fFuncValues); asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); break; @@ -1034,16 +1092,44 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2] beta = 1.0; break; case 2: // alpha != 1, beta == 1 - alpha = par[fRunInfo->GetAlphaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } beta = 1.0; break; case 3: // alpha == 1, beta != 1 alpha = 1.0; - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; case 4: // alpha != 1, beta != 1 - alpha = par[fRunInfo->GetAlphaParamNo()-1]; - beta = par[fRunInfo->GetBetaParamNo()-1]; + if (fRunInfo->GetAlphaParamNo() < MSR_PARAM_FUN_OFFSET) { // alpha is a parameter + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + } else { // alpha is function + // get function number + UInt_t funNo = fRunInfo->GetAlphaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + alpha = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } + if (fRunInfo->GetBetaParamNo() < MSR_PARAM_FUN_OFFSET) { // beta is a parameter + beta = par[fRunInfo->GetBetaParamNo()-1]; + } else { // beta is a function + // get function number + UInt_t funNo = fRunInfo->GetBetaParamNo()-MSR_PARAM_FUN_OFFSET; + // evaluate function + beta = fMsrInfo->EvalFunc(funNo, *fRunInfo->GetMap(), par); + } break; default: break;