From c215601e34e9164d14012bed6da1e64467426f8a Mon Sep 17 00:00:00 2001 From: nemu Date: Mon, 10 Mar 2008 15:49:38 +0000 Subject: [PATCH] added correct handling of alltogether unused variables --- src/ToDo.txt | 22 +- src/classes/PFitter.cpp | 12 +- src/classes/PMsrHandler.cpp | 324 ++++++++++++++++++++++ src/classes/PTheory.cpp | 22 +- src/include/PMsrHandler.h | 6 + src/tests/skewedGaussianTest/fakeData.cpp | 4 +- 6 files changed, 372 insertions(+), 18 deletions(-) diff --git a/src/ToDo.txt b/src/ToDo.txt index a7ec1d70..1b0d1e0b 100644 --- a/src/ToDo.txt +++ b/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 - -* implement table based theory functions (LF stuff) +* fix output problems of the mlog files **DONE** 08-02-22 * 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: diff --git a/src/classes/PFitter.cpp b/src/classes/PFitter.cpp index 44c3948e..713ee17d 100644 --- a/src/classes/PFitter.cpp +++ b/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; iParameterInUse(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; iParameterInUse(i) != 0)) { // 1-sigma MINOS errors std::pair err = minos(i); diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 7c592cd0..61d3bdb3 100644 --- a/src/classes/PMsrHandler.cpp +++ b/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) +//-------------------------------------------------------------------------- +/** + *

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 NOT 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) +//-------------------------------------------------------------------------- +/** + *

+ * + * \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; ifLine; + 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; iGetEntries(); i++) { + ostr = dynamic_cast(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(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; ifLine; + 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>>> 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>>> 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; posfLine; + 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(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; iGetEntries()) { + ostr = dynamic_cast(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(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; ifLine; + 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>>> 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>>> 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> fParamInUse: "; +for (unsigned int i=0; i& paramValues // val[2] = sigma-, val[3] = sigma+ - double z1 = val[3]*t/SQRT_TWO; // sigma+ - double z2 = val[2]*t/SQRT_TWO; // sigma- - double g1 = 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 zp = val[3]*t/SQRT_TWO; // sigma+ + double zm = val[2]*t/SQRT_TWO; // sigma- + double gp = TMath::Exp(-0.5*TMath::Power(t*val[3], 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 ) - - 0.5*TMath::Sin(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) * - ( - g1*2.0*z1/SQRT_PI*gsl_sf_hyperg_1F1(0.5,1.5,z1*z1) - - g2*2.0*z2/SQRT_PI*gsl_sf_hyperg_1F1(0.5,1.5,z2*z2) - ); + skg = (1.0-wp) * (gm * (TMath::Cos(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))) + + + wp * (gp * (TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*t) - + 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; } diff --git a/src/include/PMsrHandler.h b/src/include/PMsrHandler.h index 7b5a68bb..18d4e167 100644 --- a/src/include/PMsrHandler.h +++ b/src/include/PMsrHandler.h @@ -84,6 +84,8 @@ class PMsrHandler virtual double EvalFunc(unsigned int i, vector map, vector 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 ¶m); virtual bool FilterFunMapNumber(TString str, const char *filter, int &no); }; diff --git a/src/tests/skewedGaussianTest/fakeData.cpp b/src/tests/skewedGaussianTest/fakeData.cpp index 6821506f..6d5c2bff 100644 --- a/src/tests/skewedGaussianTest/fakeData.cpp +++ b/src/tests/skewedGaussianTest/fakeData.cpp @@ -51,12 +51,12 @@ typedef vector PIntVector; void fakeDataSyntax() { - cout << endl << "usage: fakeData [--dump ]"; + cout << endl << "usage: fakeData [--nemu] [--dump ]"; cout << endl << " : file holding p(B) in columns B, pB, separated by ','"; cout << endl << " comment lines start with a '#'"; cout << endl << " : holding parameters needed to generate the histograms"; cout << endl << " (see below)"; - cout << endl << " : output file name of the fake data histo file"; + cout << endl << " : output file name of the fake data histo file"; cout << endl << " --dump : option to dump p(B) and A_i(t) in a separate root file"; cout << endl; cout << endl << " structure:";