musrfit  1.9.2
PTheory.cpp
Go to the documentation of this file.
1 /***************************************************************************
2 
3  PTheory.cpp
4 
5  Author: Andreas Suter
6  e-mail: andreas.suter@psi.ch
7 
8 ***************************************************************************/
9 
10 /***************************************************************************
11  * Copyright (C) 2007-2023 by Andreas Suter *
12  * andreas.suter@psi.ch *
13  * *
14  * This program is free software; you can redistribute it and/or modify *
15  * it under the terms of the GNU General Public License as published by *
16  * the Free Software Foundation; either version 2 of the License, or *
17  * (at your option) any later version. *
18  * *
19  * This program is distributed in the hope that it will be useful, *
20  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
21  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
22  * GNU General Public License for more details. *
23  * *
24  * You should have received a copy of the GNU General Public License *
25  * along with this program; if not, write to the *
26  * Free Software Foundation, Inc., *
27  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
28  ***************************************************************************/
29 
30 #include <iostream>
31 #include <vector>
32 
33 #include <TObject.h>
34 #include <TString.h>
35 #include <TF1.h>
36 #include <TObjString.h>
37 #include <TObjArray.h>
38 #include <TClass.h>
39 #include <TMath.h>
40 
41 #include <Math/SpecFuncMathMore.h>
42 
43 #include "PMsrHandler.h"
44 #include "PTheory.h"
45 
46 #define SQRT_TWO 1.41421356237
47 #define SQRT_PI 1.77245385091
48 
49 extern std::vector<void*> gGlobalUserFcn;
50 
51 //--------------------------------------------------------------------------
52 // Constructor
53 //--------------------------------------------------------------------------
93 PTheory::PTheory(PMsrHandler *msrInfo, UInt_t runNo, const Bool_t hasParent) : fMsrInfo(msrInfo)
94 {
95  // init stuff
96  fValid = true;
97  fAdd = nullptr;
98  fMul = nullptr;
99  fUserFcnClassName = TString("");
100  fUserFcn = nullptr;
101  fDynLFdt = 0.0;
102  fSamplingTime = 0.001; // default = 1ns (units in us)
103 
104  static UInt_t lineNo = 1; // lineNo
105  static UInt_t depth = 0; // needed to handle '+' properly
106 
107  if (hasParent == false) { // reset static counters if root object
108  lineNo = 1; // the lineNo counter and the depth counter need to be
109  depth = 0; // reset for every root object (new run).
110  }
111 
112  for (UInt_t i=0; i<THEORY_MAX_PARAM; i++)
113  fPrevParam[i] = 0.0;
114 
115  // keep the number of user functions found up to this point
116  fUserFcnIdx = GetUserFcnIdx(lineNo);
117 
118  // get the input to be analyzed from the msr handler
119  PMsrLines *fullTheoryBlock = msrInfo->GetMsrTheory();
120  if (lineNo > fullTheoryBlock->size()-1) {
121  return;
122  }
123  // get the line to be parsed
124  PMsrLineStructure *line = &(*fullTheoryBlock)[lineNo];
125 
126  // copy line content to str in order to remove comments
127  TString str = line->fLine.Copy();
128 
129  // remove theory line comment if present, i.e. something starting with '('
130  Int_t index = str.Index("(");
131  if (index > 0) // theory line comment present
132  str.Resize(index);
133 
134  // remove msr-file comment if present, i.e. something starting with '#'
135  index = str.Index("#");
136  if (index > 0) // theory line comment present
137  str.Resize(index);
138 
139  // tokenize line
140  TObjArray *tokens;
141  TObjString *ostr;
142 
143  tokens = str.Tokenize(" \t");
144  if (!tokens) {
145  std::cerr << std::endl << ">> PTheory::PTheory: **SEVERE ERROR** Couldn't tokenize theory block line " << line->fLineNo << ".";
146  std::cerr << std::endl << ">> line content: " << line->fLine.Data();
147  std::cerr << std::endl;
148  exit(0);
149  }
150  ostr = dynamic_cast<TObjString*>(tokens->At(0));
151  str = ostr->GetString();
152 
153  // search the theory function
154  UInt_t idx = SearchDataBase(str);
155 
156  // function found is not defined
157  if (idx == static_cast<UInt_t>(THEORY_UNDEFINED)) {
158  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** Theory line '" << line->fLine.Data() << "'";
159  std::cerr << std::endl << ">> in line no " << line->fLineNo << " is undefined!";
160  std::cerr << std::endl;
161  fValid = false;
162  return;
163  }
164 
165  // line is a valid function, hence analyze parameters
166  if ((static_cast<UInt_t>(tokens->GetEntries()-1) < fNoOfParam) &&
167  ((idx != THEORY_USER_FCN) && (idx != THEORY_POLYNOM))) {
168  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** Theory line '" << line->fLine.Data() << "'";
169  std::cerr << std::endl << ">> in line no " << line->fLineNo;
170  std::cerr << std::endl << ">> expecting " << fgTheoDataBase[idx].fNoOfParam << ", but found " << tokens->GetEntries()-1;
171  std::cerr << std::endl;
172  fValid = false;
173  }
174  // keep function index
175  fType = idx;
176  // filter out the parameters
177  Int_t status;
178  UInt_t value;
179  Bool_t ok = false;
180  for (Int_t i=1; i<tokens->GetEntries(); i++) {
181  ostr = dynamic_cast<TObjString*>(tokens->At(i));
182  str = ostr->GetString();
183 
184  // if userFcn, the first entry is the function name and needs to be handled specially
185  if ((fType == THEORY_USER_FCN) && ((i == 1) || (i == 2))) {
186  if (i == 1) {
187  fUserFcnSharedLibName = str;
188  }
189  if (i == 2) {
190  fUserFcnClassName = str;
191  }
192  continue;
193  }
194 
195  // check if str is map
196  if (str.Contains("map")) {
197  status = sscanf(str.Data(), "map%u", &value);
198  if (status == 1) { // everthing ok
199  ok = true;
200  // get parameter from map
201  PIntVector maps = *(*msrInfo->GetMsrRunList())[runNo].GetMap();
202  if ((value <= maps.size()) && (value > 0)) { // everything fine
203  fParamNo.push_back(maps[value-1]-1);
204  } else { // map index out of range
205  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** map index " << value << " out of range! See line no " << line->fLineNo;
206  std::cerr << std::endl;
207  fValid = false;
208  }
209  } else { // something wrong
210  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR**: map '" << str.Data() << "' not allowed. See line no " << line->fLineNo;
211  std::cerr << std::endl;
212  fValid = false;
213  }
214  } else if (str.Contains("fun")) { // check if str is fun
215  status = sscanf(str.Data(), "fun%u", &value);
216  if (status == 1) { // everthing ok
217  ok = true;
218  // handle function, i.e. get, from the function number x (FUNx), the function index,
219  // add function offset and fill "parameter vector"
220  fParamNo.push_back(msrInfo->GetFuncIndex(value)+MSR_PARAM_FUN_OFFSET);
221  } else { // something wrong
222  fValid = false;
223  }
224  } else { // check if str is param no
225  status = sscanf(str.Data(), "%u", &value);
226  if (status == 1) { // everthing ok
227  ok = true;
228  fParamNo.push_back(value-1);
229  }
230  // check if one of the valid entries was found
231  if (!ok) {
232  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** '" << str.Data() << "' not allowed. See line no " << line->fLineNo;
233  std::cerr << std::endl;
234  fValid = false;
235  }
236  }
237  }
238 
239  // call the next line (only if valid so far and not the last line)
240  // handle '*'
241  if (fValid && (lineNo < fullTheoryBlock->size()-1)) {
242  line = &(*fullTheoryBlock)[lineNo+1];
243  if (!line->fLine.Contains("+")) { // make sure next line is not a '+'
244  depth++;
245  lineNo++;
246  fMul = new PTheory(msrInfo, runNo, true);
247  depth--;
248  }
249  }
250  // call the next line (only if valid so far and not the last line)
251  // handle '+'
252  if (fValid && (lineNo < fullTheoryBlock->size()-1)) {
253  line = &(*fullTheoryBlock)[lineNo+1];
254  if ((depth == 0) && line->fLine.Contains("+")) {
255  lineNo += 2; // go to the next theory function line
256  fAdd = new PTheory(msrInfo, runNo, true);
257  }
258  }
259 
260  // make clean and tidy theory block for the msr-file
261  if (fValid && !hasParent) { // parent theory object
262  MakeCleanAndTidyTheoryBlock(fullTheoryBlock);
263  }
264 
265  // check if user function, if so, check if it is reachable (root) and if yes invoke object
266  if (!fUserFcnClassName.IsWhitespace()) {
267  std::cout << std::endl << ">> user function class name: " << fUserFcnClassName.Data() << std::endl;
268  if (!TClass::GetDict(fUserFcnClassName.Data())) {
269  if (gSystem->Load(fUserFcnSharedLibName.Data()) < 0) {
270  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** user function class '" << fUserFcnClassName.Data() << "' not found.";
271  std::cerr << std::endl << ">> Tried to load " << fUserFcnSharedLibName.Data() << " but failed.";
272  std::cerr << std::endl << ">> See line no " << line->fLineNo;
273  std::cerr << std::endl;
274  fValid = false;
275  // clean up
276  if (tokens) {
277  delete tokens;
278  tokens = nullptr;
279  }
280  return;
281  } else if (!TClass::GetDict(fUserFcnClassName.Data())) {
282  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** user function class '" << fUserFcnClassName.Data() << "' not found.";
283  std::cerr << std::endl << ">> " << fUserFcnSharedLibName.Data() << " loaded successfully, but no dictionary present.";
284  std::cerr << std::endl << ">> See line no " << line->fLineNo;
285  std::cerr << std::endl;
286  fValid = false;
287  // clean up
288  if (tokens) {
289  delete tokens;
290  tokens = nullptr;
291  }
292  return;
293  }
294  }
295 
296  // invoke user function object
297  fUserFcn = nullptr;
298  fUserFcn = static_cast<PUserFcnBase*>(TClass::GetClass(fUserFcnClassName.Data())->New());
299  if (fUserFcn == nullptr) {
300  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** user function object could not be invoked. See line no " << line->fLineNo;
301  std::cerr << std::endl;
302  fValid = false;
303  return;
304  } else { // user function valid, hence expand the fUserParam vector to the proper size
305  fUserParam.resize(fParamNo.size());
306  }
307 
308  // check if the global part of the user function is needed
309  if (fUserFcn->NeedGlobalPart()) {
310  fUserFcn->SetGlobalPart(gGlobalUserFcn, fUserFcnIdx); // invoke or retrieve global user function object
311  if (!fUserFcn->GlobalPartIsValid()) {
312  std::cerr << std::endl << ">> PTheory::PTheory: **ERROR** global user function object could not be invoked/retrived. See line no " << line->fLineNo;
313  std::cerr << std::endl;
314  fValid = false;
315  }
316  }
317  }
318 
319  // clean up
320  if (tokens) {
321  delete tokens;
322  tokens = nullptr;
323  }
324 }
325 
326 //--------------------------------------------------------------------------
327 // Destructor
328 //--------------------------------------------------------------------------
332 PTheory::~PTheory()
333 {
334  fParamNo.clear();
335  fUserParam.clear();
336 
337  fLFIntegral.clear();
338  fDynLFFuncValue.clear();
339 
340  // recursive clean up
341  CleanUp(this);
342 
343  if (fUserFcn) {
344  delete fUserFcn;
345  fUserFcn = nullptr;
346  }
347 
348  gGlobalUserFcn.clear();
349 }
350 
351 //--------------------------------------------------------------------------
352 // IsValid
353 //--------------------------------------------------------------------------
361 Bool_t PTheory::IsValid()
362 {
363 
364  if (fMul) {
365  if (fAdd) {
366  return (fValid && fMul->IsValid() && fAdd->IsValid());
367  } else {
368  return (fValid && fMul->IsValid());
369  }
370  } else {
371  if (fAdd) {
372  return (fValid && fAdd->IsValid());
373  } else {
374  return fValid;
375  }
376  }
377 }
378 
379 //--------------------------------------------------------------------------
390 Double_t PTheory::Func(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
391 {
392  if (fMul) {
393  if (fAdd) { // fMul != 0 && fAdd != 0
394  switch (fType) {
395  case THEORY_CONST:
396  return Constant(paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
397  fAdd->Func(t, paramValues, funcValues);
398  case THEORY_ASYMMETRY:
399  return Asymmetry(paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
400  fAdd->Func(t, paramValues, funcValues);
401  case THEORY_SIMPLE_EXP:
402  return SimpleExp(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
403  fAdd->Func(t, paramValues, funcValues);
404  case THEORY_GENERAL_EXP:
405  return GeneralExp(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
406  fAdd->Func(t, paramValues, funcValues);
407  case THEORY_SIMPLE_GAUSS:
408  return SimpleGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
409  fAdd->Func(t, paramValues, funcValues);
410  case THEORY_STATIC_GAUSS_KT:
411  return StaticGaussKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
412  fAdd->Func(t, paramValues, funcValues);
413  case THEORY_STATIC_GAUSS_KT_LF:
414  return StaticGaussKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
415  fAdd->Func(t, paramValues, funcValues);
416  case THEORY_DYNAMIC_GAUSS_KT_LF:
417  return DynamicGaussKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
418  fAdd->Func(t, paramValues, funcValues);
419  case THEORY_STATIC_LORENTZ_KT:
420  return StaticLorentzKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
421  fAdd->Func(t, paramValues, funcValues);
422  case THEORY_STATIC_LORENTZ_KT_LF:
423  return StaticLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
424  fAdd->Func(t, paramValues, funcValues);
425  case THEORY_DYNAMIC_LORENTZ_KT_LF:
426  return DynamicLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
427  fAdd->Func(t, paramValues, funcValues);
428  case THEORY_COMBI_LGKT:
429  return CombiLGKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
430  fAdd->Func(t, paramValues, funcValues);
431  case THEORY_STR_KT:
432  return StrKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
433  fAdd->Func(t, paramValues, funcValues);
434  case THEORY_SPIN_GLASS:
435  return SpinGlass(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
436  fAdd->Func(t, paramValues, funcValues);
437  case THEORY_RANDOM_ANISOTROPIC_HYPERFINE:
438  return RandomAnisotropicHyperfine(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
439  fAdd->Func(t, paramValues, funcValues);
440  case THEORY_ABRAGAM:
441  return Abragam(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
442  fAdd->Func(t, paramValues, funcValues);
443  case THEORY_TF_COS:
444  return TFCos(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
445  fAdd->Func(t, paramValues, funcValues);
446  case THEORY_INTERNAL_FIELD:
447  return InternalField(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
448  fAdd->Func(t, paramValues, funcValues);
449  case THEORY_INTERNAL_FIELD_KORNILOV:
450  return InternalFieldGK(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
451  fAdd->Func(t, paramValues, funcValues);
452  case THEORY_INTERNAL_FIELD_LARKIN:
453  return InternalFieldLL(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
454  fAdd->Func(t, paramValues, funcValues);
455  case THEORY_BESSEL:
456  return Bessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
457  fAdd->Func(t, paramValues, funcValues);
458  case THEORY_INTERNAL_BESSEL:
459  return InternalBessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
460  fAdd->Func(t, paramValues, funcValues);
461  case THEORY_SKEWED_GAUSS:
462  return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
463  fAdd->Func(t, paramValues, funcValues);
464  case THEORY_STATIC_ZF_NK:
465  return StaticNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
466  fAdd->Func(t, paramValues, funcValues);
467  case THEORY_STATIC_TF_NK:
468  return StaticNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
469  fAdd->Func(t, paramValues, funcValues);
470  case THEORY_DYNAMIC_ZF_NK:
471  return DynamicNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
472  fAdd->Func(t, paramValues, funcValues);
473  case THEORY_DYNAMIC_TF_NK:
474  return DynamicNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
475  fAdd->Func(t, paramValues, funcValues);
476  case THEORY_POLYNOM:
477  return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
478  fAdd->Func(t, paramValues, funcValues);
479  case THEORY_MU_MINUS_EXP:
480  return MuMinusExpTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
481  fAdd->Func(t, paramValues, funcValues);
482  case THEORY_USER_FCN:
483  return UserFcn(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues) +
484  fAdd->Func(t, paramValues, funcValues);
485  default:
486  std::cerr << std::endl << ">> PTheory::Func: **PANIC ERROR** You never should have reached this line?!?! (" << fType << ")";
487  std::cerr << std::endl;
488  exit(0);
489  }
490  } else { // fMul !=0 && fAdd == 0
491  switch (fType) {
492  case THEORY_CONST:
493  return Constant(paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
494  case THEORY_ASYMMETRY:
495  return Asymmetry(paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
496  case THEORY_SIMPLE_EXP:
497  return SimpleExp(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
498  case THEORY_GENERAL_EXP:
499  return GeneralExp(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
500  case THEORY_SIMPLE_GAUSS:
501  return SimpleGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
502  case THEORY_STATIC_GAUSS_KT:
503  return StaticGaussKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
504  case THEORY_STATIC_GAUSS_KT_LF:
505  return StaticGaussKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
506  case THEORY_DYNAMIC_GAUSS_KT_LF:
507  return DynamicGaussKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
508  case THEORY_STATIC_LORENTZ_KT:
509  return StaticLorentzKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
510  case THEORY_STATIC_LORENTZ_KT_LF:
511  return StaticLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
512  case THEORY_DYNAMIC_LORENTZ_KT_LF:
513  return DynamicLorentzKTLF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
514  case THEORY_COMBI_LGKT:
515  return CombiLGKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
516  case THEORY_STR_KT:
517  return StrKT(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
518  case THEORY_SPIN_GLASS:
519  return SpinGlass(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
520  case THEORY_RANDOM_ANISOTROPIC_HYPERFINE:
521  return RandomAnisotropicHyperfine(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
522  case THEORY_ABRAGAM:
523  return Abragam(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
524  case THEORY_TF_COS:
525  return TFCos(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
526  case THEORY_INTERNAL_FIELD:
527  return InternalField(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
528  case THEORY_INTERNAL_FIELD_KORNILOV:
529  return InternalFieldGK(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
530  case THEORY_INTERNAL_FIELD_LARKIN:
531  return InternalFieldLL(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
532  case THEORY_BESSEL:
533  return Bessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
534  case THEORY_INTERNAL_BESSEL:
535  return InternalBessel(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
536  case THEORY_SKEWED_GAUSS:
537  return SkewedGauss(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
538  case THEORY_STATIC_ZF_NK:
539  return StaticNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
540  case THEORY_STATIC_TF_NK:
541  return StaticNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
542  case THEORY_DYNAMIC_ZF_NK:
543  return DynamicNKZF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
544  case THEORY_DYNAMIC_TF_NK:
545  return DynamicNKTF (t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
546  case THEORY_MU_MINUS_EXP:
547  return MuMinusExpTF(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
548  case THEORY_POLYNOM:
549  return Polynom(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
550  case THEORY_USER_FCN:
551  return UserFcn(t, paramValues, funcValues) * fMul->Func(t, paramValues, funcValues);
552  default:
553  std::cerr << std::endl << ">> PTheory::Func: **PANIC ERROR** You never should have reached this line?!?! (" << fType << ")";
554  std::cerr << std::endl;
555  exit(0);
556  }
557  }
558  } else { // fMul == 0 && fAdd != 0
559  if (fAdd) {
560  switch (fType) {
561  case THEORY_CONST:
562  return Constant(paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
563  case THEORY_ASYMMETRY:
564  return Asymmetry(paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
565  case THEORY_SIMPLE_EXP:
566  return SimpleExp(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
567  case THEORY_GENERAL_EXP:
568  return GeneralExp(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
569  case THEORY_SIMPLE_GAUSS:
570  return SimpleGauss(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
571  case THEORY_STATIC_GAUSS_KT:
572  return StaticGaussKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
573  case THEORY_STATIC_GAUSS_KT_LF:
574  return StaticGaussKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
575  case THEORY_DYNAMIC_GAUSS_KT_LF:
576  return DynamicGaussKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
577  case THEORY_STATIC_LORENTZ_KT:
578  return StaticLorentzKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
579  case THEORY_STATIC_LORENTZ_KT_LF:
580  return StaticLorentzKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
581  case THEORY_DYNAMIC_LORENTZ_KT_LF:
582  return DynamicLorentzKTLF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
583  case THEORY_COMBI_LGKT:
584  return CombiLGKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
585  case THEORY_STR_KT:
586  return StrKT(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
587  case THEORY_SPIN_GLASS:
588  return SpinGlass(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
589  case THEORY_RANDOM_ANISOTROPIC_HYPERFINE:
590  return RandomAnisotropicHyperfine(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
591  case THEORY_ABRAGAM:
592  return Abragam(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
593  case THEORY_TF_COS:
594  return TFCos(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
595  case THEORY_INTERNAL_FIELD:
596  return InternalField(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
597  case THEORY_INTERNAL_FIELD_KORNILOV:
598  return InternalFieldGK(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
599  case THEORY_INTERNAL_FIELD_LARKIN:
600  return InternalFieldLL(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
601  case THEORY_BESSEL:
602  return Bessel(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
603  case THEORY_INTERNAL_BESSEL:
604  return InternalBessel(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
605  case THEORY_SKEWED_GAUSS:
606  return SkewedGauss(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
607  case THEORY_STATIC_ZF_NK:
608  return StaticNKZF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
609  case THEORY_STATIC_TF_NK:
610  return StaticNKTF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
611  case THEORY_DYNAMIC_ZF_NK:
612  return DynamicNKZF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
613  case THEORY_DYNAMIC_TF_NK:
614  return DynamicNKTF (t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
615  case THEORY_MU_MINUS_EXP:
616  return MuMinusExpTF(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
617  case THEORY_POLYNOM:
618  return Polynom(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
619  case THEORY_USER_FCN:
620  return UserFcn(t, paramValues, funcValues) + fAdd->Func(t, paramValues, funcValues);
621  default:
622  std::cerr << std::endl << ">> PTheory::Func: **PANIC ERROR** You never should have reached this line?!?! (" << fType << ")";
623  std::cerr << std::endl;
624  exit(0);
625  }
626  } else { // fMul == 0 && fAdd == 0
627  switch (fType) {
628  case THEORY_CONST:
629  return Constant(paramValues, funcValues);
630  case THEORY_ASYMMETRY:
631  return Asymmetry(paramValues, funcValues);
632  case THEORY_SIMPLE_EXP:
633  return SimpleExp(t, paramValues, funcValues);
634  case THEORY_GENERAL_EXP:
635  return GeneralExp(t, paramValues, funcValues);
636  case THEORY_SIMPLE_GAUSS:
637  return SimpleGauss(t, paramValues, funcValues);
638  case THEORY_STATIC_GAUSS_KT:
639  return StaticGaussKT(t, paramValues, funcValues);
640  case THEORY_STATIC_GAUSS_KT_LF:
641  return StaticGaussKTLF(t, paramValues, funcValues);
642  case THEORY_DYNAMIC_GAUSS_KT_LF:
643  return DynamicGaussKTLF(t, paramValues, funcValues);
644  case THEORY_STATIC_LORENTZ_KT:
645  return StaticLorentzKT(t, paramValues, funcValues);
646  case THEORY_STATIC_LORENTZ_KT_LF:
647  return StaticLorentzKTLF(t, paramValues, funcValues);
648  case THEORY_DYNAMIC_LORENTZ_KT_LF:
649  return DynamicLorentzKTLF(t, paramValues, funcValues);
650  case THEORY_COMBI_LGKT:
651  return CombiLGKT(t, paramValues, funcValues);
652  case THEORY_STR_KT:
653  return StrKT(t, paramValues, funcValues);
654  case THEORY_SPIN_GLASS:
655  return SpinGlass(t, paramValues, funcValues);
656  case THEORY_RANDOM_ANISOTROPIC_HYPERFINE:
657  return RandomAnisotropicHyperfine(t, paramValues, funcValues);
658  case THEORY_ABRAGAM:
659  return Abragam(t, paramValues, funcValues);
660  case THEORY_TF_COS:
661  return TFCos(t, paramValues, funcValues);
662  case THEORY_INTERNAL_FIELD:
663  return InternalField(t, paramValues, funcValues);
664  case THEORY_INTERNAL_FIELD_KORNILOV:
665  return InternalFieldGK(t, paramValues, funcValues);
666  case THEORY_INTERNAL_FIELD_LARKIN:
667  return InternalFieldLL(t, paramValues, funcValues);
668  case THEORY_BESSEL:
669  return Bessel(t, paramValues, funcValues);
670  case THEORY_INTERNAL_BESSEL:
671  return InternalBessel(t, paramValues, funcValues);
672  case THEORY_SKEWED_GAUSS:
673  return SkewedGauss(t, paramValues, funcValues);
674  case THEORY_STATIC_ZF_NK:
675  return StaticNKZF(t, paramValues, funcValues);
676  case THEORY_STATIC_TF_NK:
677  return StaticNKTF(t, paramValues, funcValues);
678  case THEORY_DYNAMIC_ZF_NK:
679  return DynamicNKZF(t, paramValues, funcValues);
680  case THEORY_DYNAMIC_TF_NK:
681  return DynamicNKTF(t, paramValues, funcValues);
682  case THEORY_MU_MINUS_EXP:
683  return MuMinusExpTF(t, paramValues, funcValues);
684  case THEORY_POLYNOM:
685  return Polynom(t, paramValues, funcValues);
686  case THEORY_USER_FCN:
687  return UserFcn(t, paramValues, funcValues);
688  default:
689  std::cerr << std::endl << ">> PTheory::Func: **PANIC ERROR** You never should have reached this line?!?! (" << fType << ")";
690  std::cerr << std::endl;
691  exit(0);
692  }
693  }
694  }
695 }
696 
697 //--------------------------------------------------------------------------
716 void PTheory::CleanUp(PTheory *theo)
717 {
718  if (theo->fMul) { // '*' present
719  delete theo->fMul;
720  theo->fMul = nullptr;
721  }
722 
723  if (theo->fAdd) {
724  delete theo->fAdd;
725  theo->fAdd = nullptr;
726  }
727 }
728 
729 //--------------------------------------------------------------------------
739 Int_t PTheory::SearchDataBase(TString name)
740 {
741  Int_t idx = THEORY_UNDEFINED;
742 
743  for (UInt_t i=0; i<THEORY_MAX; i++) {
744  if (!name.CompareTo(fgTheoDataBase[i].fName, TString::kIgnoreCase) ||
745  !name.CompareTo(fgTheoDataBase[i].fAbbrev, TString::kIgnoreCase)) {
746  idx = static_cast<Int_t>(fgTheoDataBase[i].fType);
747  fType = fgTheoDataBase[i].fType;
748  fNoOfParam = fgTheoDataBase[i].fNoOfParam;
749  }
750  }
751 
752  return idx;
753 }
754 
755 //--------------------------------------------------------------------------
756 // GetUserFcnIdx (private)
757 //--------------------------------------------------------------------------
765 Int_t PTheory::GetUserFcnIdx(UInt_t lineNo) const
766 {
767  Int_t userFcnIdx = -1;
768 
769  // retrieve the theory block from the msr-file handler
770  PMsrLines *fullTheoryBlock = fMsrInfo->GetMsrTheory();
771 
772  // make sure that lineNo is within proper bounds
773  if (lineNo > fullTheoryBlock->size())
774  return userFcnIdx;
775 
776  // count the number of user function present up to the lineNo
777  for (UInt_t i=1; i<=lineNo; i++) {
778  if (fullTheoryBlock->at(i).fLine.Contains("userFcn", TString::kIgnoreCase))
779  userFcnIdx++;
780  }
781 
782  return userFcnIdx;
783 }
784 
785 //--------------------------------------------------------------------------
786 // MakeCleanAndTidyTheoryBlock private
787 //--------------------------------------------------------------------------
794 void PTheory::MakeCleanAndTidyTheoryBlock(PMsrLines *fullTheoryBlock)
795 {
796  PMsrLineStructure *line;
797  TString str, tidy;
798  Char_t substr[256];
799  TObjArray *tokens = nullptr;
800  TObjString *ostr = nullptr;
801  Int_t idx = THEORY_UNDEFINED;
802 
803  for (UInt_t i=1; i<fullTheoryBlock->size(); i++) {
804  // get the line to be prettyfied
805  line = &(*fullTheoryBlock)[i];
806  // copy line content to str in order to remove comments
807  str = line->fLine.Copy();
808  // remove theory line comment if present, i.e. something starting with '('
809  Int_t index = str.Index("(");
810  if (index > 0) // theory line comment present
811  str.Resize(index);
812  // tokenize line
813  tokens = str.Tokenize(" \t");
814  // make a handable string out of the asymmetry token
815  ostr = dynamic_cast<TObjString*>(tokens->At(0));
816  str = ostr->GetString();
817  // check if the line is just a '+' if so nothing to be done
818  if (str.Contains("+"))
819  continue;
820  // check if the function is a polynom
821  if (!str.CompareTo("p") || str.Contains("polynom")) {
822  MakeCleanAndTidyPolynom(i, fullTheoryBlock);
823  continue;
824  }
825  // check if the function is a userFcn
826  if (!str.CompareTo("u") || str.Contains("userFcn")) {
827  MakeCleanAndTidyUserFcn(i, fullTheoryBlock);
828  continue;
829  }
830  // search the theory function
831  for (UInt_t j=0; j<THEORY_MAX; j++) {
832  if (!str.CompareTo(fgTheoDataBase[j].fName, TString::kIgnoreCase) ||
833  !str.CompareTo(fgTheoDataBase[j].fAbbrev, TString::kIgnoreCase)) {
834  idx = static_cast<Int_t>(fgTheoDataBase[j].fType);
835  }
836  }
837  // check if theory is indeed defined. This should not be necessay at this point but ...
838  if (idx == THEORY_UNDEFINED)
839  return;
840  // check that there enough tokens. This should not be necessay at this point but ...
841  if (static_cast<UInt_t>(tokens->GetEntries()) < fgTheoDataBase[idx].fNoOfParam + 1)
842  return;
843  // make tidy string
844  snprintf(substr, sizeof(substr), "%-10s", fgTheoDataBase[idx].fName.Data());
845  tidy = TString(substr);
846  for (Int_t j=1; j<tokens->GetEntries(); j++) {
847  ostr = dynamic_cast<TObjString*>(tokens->At(j));
848  str = ostr->GetString();
849  snprintf(substr, sizeof(substr), "%6s", str.Data());
850  tidy += TString(substr);
851  }
852  if (fgTheoDataBase[idx].fComment.Length() != 0) {
853  if (tidy.Length() < 35) {
854  for (Int_t k=0; k<35-tidy.Length(); k++)
855  tidy += TString(" ");
856  } else {
857  tidy += TString(" ");
858  }
859  if (static_cast<UInt_t>(tokens->GetEntries()) == fgTheoDataBase[idx].fNoOfParam + 1) // no tshift
860  tidy += fgTheoDataBase[idx].fComment;
861  else
862  tidy += fgTheoDataBase[idx].fCommentTimeShift;
863  }
864  // write tidy string back into theory block
865  (*fullTheoryBlock)[i].fLine = tidy;
866 
867  // clean up
868  if (tokens) {
869  delete tokens;
870  tokens = nullptr;
871  }
872  }
873 
874 }
875 
876 //--------------------------------------------------------------------------
877 // MakeCleanAndTidyPolynom private
878 //--------------------------------------------------------------------------
885 void PTheory::MakeCleanAndTidyPolynom(UInt_t i, PMsrLines *fullTheoryBlock)
886 {
887  PMsrLineStructure *line;
888  TString str, tidy;
889  TObjArray *tokens = nullptr;
890  TObjString *ostr;
891  Char_t substr[256];
892 
893  // init tidy
894  tidy = TString("polynom ");
895  // get the line to be prettyfied
896  line = &(*fullTheoryBlock)[i];
897  // copy line content to str in order to remove comments
898  str = line->fLine.Copy();
899  // tokenize line
900  tokens = str.Tokenize(" \t");
901 
902  // check if comment is already present, and if yes ignore it by setting max correctly
903  Int_t max = tokens->GetEntries();
904  for (Int_t j=1; j<max; j++) {
905  ostr = dynamic_cast<TObjString*>(tokens->At(j));
906  str = ostr->GetString();
907  if (str.Contains("(")) { // comment present
908  max=j;
909  break;
910  }
911  }
912 
913  for (Int_t j=1; j<max; j++) {
914  ostr = dynamic_cast<TObjString*>(tokens->At(j));
915  str = ostr->GetString();
916  snprintf(substr, sizeof(substr), "%6s", str.Data());
917  tidy += TString(substr);
918  }
919 
920  // add comment
921  tidy += " (tshift p0 p1 ... pn)";
922 
923  // write tidy string back into theory block
924  (*fullTheoryBlock)[i].fLine = tidy;
925 
926  // clean up
927  if (tokens) {
928  delete tokens;
929  tokens = nullptr;
930  }
931 }
932 
933 //--------------------------------------------------------------------------
934 // MakeCleanAndTidyUserFcn private
935 //--------------------------------------------------------------------------
942 void PTheory::MakeCleanAndTidyUserFcn(UInt_t i, PMsrLines *fullTheoryBlock)
943 {
944  PMsrLineStructure *line;
945  TString str, tidy;
946  TObjArray *tokens = nullptr;
947  TObjString *ostr;
948 
949  // init tidy
950  tidy = TString("userFcn ");
951  // get the line to be prettyfied
952  line = &(*fullTheoryBlock)[i];
953  // copy line content to str in order to remove comments
954  str = line->fLine.Copy();
955  // tokenize line
956  tokens = str.Tokenize(" \t");
957 
958  for (Int_t j=1; j<tokens->GetEntries(); j++) {
959  ostr = dynamic_cast<TObjString*>(tokens->At(j));
960  str = ostr->GetString();
961  tidy += TString(" ") + str;
962  }
963 
964  // write tidy string back into theory block
965  (*fullTheoryBlock)[i].fLine = tidy;
966 
967  // clean up
968  if (tokens) {
969  delete tokens;
970  tokens = nullptr;
971  }
972 }
973 
974 //--------------------------------------------------------------------------
986 Double_t PTheory::Constant(const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
987 {
988  // expected parameters: const
989 
990  Double_t constant;
991 
992  // check if FUNCTIONS are used
993  if (fParamNo[0] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
994  constant = paramValues[fParamNo[0]];
995  } else {
996  constant = funcValues[fParamNo[0]-MSR_PARAM_FUN_OFFSET];
997  }
998 
999  return constant;
1000 }
1001 
1002 //--------------------------------------------------------------------------
1014 Double_t PTheory::Asymmetry(const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1015 {
1016  // expected parameters: asym
1017 
1018  Double_t asym;
1019 
1020  // check if FUNCTIONS are used
1021  if (fParamNo[0] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1022  asym = paramValues[fParamNo[0]];
1023  } else {
1024  asym = funcValues[fParamNo[0]-MSR_PARAM_FUN_OFFSET];
1025  }
1026 
1027  return asym;
1028 }
1029 
1030 //--------------------------------------------------------------------------
1044 Double_t PTheory::SimpleExp(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1045 {
1046  // expected parameters: lambda [tshift]
1047 
1048  Double_t val[2];
1049 
1050  assert(fParamNo.size() <= 2);
1051 
1052  // check if FUNCTIONS are used
1053  for (UInt_t i=0; i<fParamNo.size(); i++) {
1054  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1055  val[i] = paramValues[fParamNo[i]];
1056  } else { // function
1057  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1058  }
1059  }
1060 
1061  Double_t tt;
1062  if (fParamNo.size() == 1) // no tshift
1063  tt = t;
1064  else // tshift present
1065  tt = t-val[1];
1066 
1067  return TMath::Exp(-tt*val[0]);
1068 }
1069 
1070 //--------------------------------------------------------------------------
1084 Double_t PTheory::GeneralExp(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1085 {
1086  // expected parameters: lambda beta [tshift]
1087 
1088  Double_t val[3];
1089  Double_t result;
1090 
1091  assert(fParamNo.size() <= 3);
1092 
1093  // check if FUNCTIONS are used
1094  for (UInt_t i=0; i<fParamNo.size(); i++) {
1095  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1096  val[i] = paramValues[fParamNo[i]];
1097  } else { // function
1098  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1099  }
1100  }
1101 
1102  Double_t tt;
1103  if (fParamNo.size() == 2) // no tshift
1104  tt = t;
1105  else // tshift present
1106  tt = t-val[2];
1107 
1108  // check if tt*val[0] < 0 and
1109  if ((tt*val[0] < 0) && (trunc(val[1])-val[1] != 0.0)) {
1110  result = 0.0;
1111  } else {
1112  result = TMath::Exp(-TMath::Power(tt*val[0], val[1]));
1113  }
1114 
1115  return result;
1116 }
1117 
1118 //--------------------------------------------------------------------------
1132 Double_t PTheory::SimpleGauss(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1133 {
1134  // expected parameters: sigma [tshift]
1135 
1136  Double_t val[2];
1137 
1138  assert(fParamNo.size() <= 2);
1139 
1140  // check if FUNCTIONS are used
1141  for (UInt_t i=0; i<fParamNo.size(); i++) {
1142  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1143  val[i] = paramValues[fParamNo[i]];
1144  } else { // function
1145  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1146  }
1147  }
1148 
1149  Double_t tt;
1150  if (fParamNo.size() == 1) // no tshift
1151  tt = t;
1152  else // tshift present
1153  tt = t-val[1];
1154 
1155  return TMath::Exp(-0.5*TMath::Power(tt*val[0], 2.0));
1156 }
1157 
1158 //--------------------------------------------------------------------------
1172 Double_t PTheory::StaticGaussKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1173 {
1174  // expected parameters: sigma [tshift]
1175 
1176  Double_t val[2];
1177 
1178  assert(fParamNo.size() <= 2);
1179 
1180  // check if FUNCTIONS are used
1181  for (UInt_t i=0; i<fParamNo.size(); i++) {
1182  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1183  val[i] = paramValues[fParamNo[i]];
1184  } else { // function
1185  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1186  }
1187  }
1188 
1189  Double_t sigma_t_2;
1190  if (fParamNo.size() == 1) // no tshift
1191  sigma_t_2 = t*t*val[0]*val[0];
1192  else // tshift present
1193  sigma_t_2 = (t-val[1])*(t-val[1])*val[0]*val[0];
1194 
1195  return 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
1196 }
1197 
1198 //--------------------------------------------------------------------------
1214 Double_t PTheory::StaticGaussKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1215 {
1216 
1217  // expected parameters: frequency damping [tshift]
1218 
1219  Double_t val[3];
1220  Double_t result;
1221 
1222  assert(fParamNo.size() <= 3);
1223 
1224  // check if FUNCTIONS are used
1225  for (UInt_t i=0; i<fParamNo.size(); i++) {
1226  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1227  val[i] = paramValues[fParamNo[i]];
1228  } else { // function
1229  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1230  }
1231  }
1232 
1233  // check if all parameters == 0
1234  if ((val[0] == 0.0) && (val[1] == 0.0))
1235  return 1.0;
1236 
1237  // check if the parameter values have changed, and if yes recalculate the non-analytic integral
1238  // check only the first two parameters since the tshift is irrelevant for the LF-integral calculation!!
1239  Bool_t newParam = false;
1240  for (UInt_t i=0; i<2; i++) {
1241  if (val[i] != fPrevParam[i]) {
1242  newParam = true;
1243  break;
1244  }
1245  }
1246 
1247  if (newParam) { // new parameters found
1248  {
1249  for (UInt_t i=0; i<2; i++)
1250  fPrevParam[i] = val[i];
1251  CalculateGaussLFIntegral(val);
1252  }
1253  }
1254 
1255  Double_t tt;
1256  if (fParamNo.size() == 2) // no tshift
1257  tt = t;
1258  else // tshift present
1259  tt = t-val[2];
1260 
1261  if (tt < 0.0) // for times < 0 return a function value of 1.0
1262  return 1.0;
1263 
1264  Double_t sigma_t_2 = 0.0;
1265  if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use the ZF formula
1266  sigma_t_2 = tt*tt*val[1]*val[1];
1267  result = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
1268  } else if (val[1]/val[0] > 79.5775) { // check if Delta/w0 > 500.0, in which case the ZF formula is used
1269  sigma_t_2 = tt*tt*val[1]*val[1];
1270  result = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
1271  } else {
1272  Double_t delta = val[1];
1273  Double_t w0 = 2.0*TMath::Pi()*val[0];
1274 
1275  result = 1.0 - 2.0*TMath::Power(delta/w0,2.0)*(1.0 -
1276  TMath::Exp(-0.5*TMath::Power(delta*tt, 2.0))*TMath::Cos(w0*tt)) +
1277  GetLFIntegralValue(tt);
1278  }
1279 
1280  return result;
1281 
1282 }
1283 
1284 //--------------------------------------------------------------------------
1303 Double_t PTheory::DynamicGaussKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1304 {
1305  // expected parameters: frequency damping hopping [tshift]
1306 
1307  Double_t val[4];
1308  Double_t result = 0.0;
1309  Bool_t useKeren = false;
1310 
1311  assert(fParamNo.size() <= 4);
1312 
1313  // check if FUNCTIONS are used
1314  for (UInt_t i=0; i<fParamNo.size(); i++) {
1315  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1316  val[i] = paramValues[fParamNo[i]];
1317  } else { // function
1318  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1319  }
1320  }
1321 
1322  // check if all parameters == 0
1323  if ((val[0] == 0.0) && (val[1] == 0.0) && (val[2] == 0.0))
1324  return 1.0;
1325 
1326  // make sure that damping and hopping are positive definite
1327  if (val[1] < 0.0)
1328  val[1] = -val[1];
1329  if (val[2] < 0.0)
1330  val[2] = -val[2];
1331 
1332  // check that Delta != 0, if not (i.e. stupid parameter) return 1, which is the correct limit
1333  if (fabs(val[1]) < 1.0e-6) {
1334  return 1.0;
1335  }
1336 
1337  // check if Keren approximation can be used
1338  if (val[2]/val[1] > 5.0) // nu/Delta > 5.0
1339  useKeren = true;
1340 
1341  if (!useKeren) {
1342  // check if the parameter values have changed, and if yes recalculate the non-analytic integral
1343  // check only the first three parameters since the tshift is irrelevant for the LF-integral calculation!!
1344  Bool_t newParam = false;
1345  for (UInt_t i=0; i<3; i++) {
1346  if (val[i] != fPrevParam[i]) {
1347  newParam = true;
1348  break;
1349  }
1350  }
1351 
1352  if (newParam) { // new parameters found
1353  for (UInt_t i=0; i<3; i++)
1354  fPrevParam[i] = val[i];
1355  CalculateDynKTLF(val, 0); // 0 means Gauss
1356  }
1357  }
1358 
1359  Double_t tt;
1360  if (fParamNo.size() == 3) // no tshift
1361  tt = t;
1362  else // tshift present
1363  tt = t-val[3];
1364 
1365  if (tt < 0.0) // for times < 0 return a function value of 0.0
1366  return 0.0;
1367 
1368 
1369  if (useKeren) { // see PRB50, 10039 (1994)
1370  Double_t wL = TWO_PI * val[0];
1371  Double_t wL2 = wL*wL;
1372  Double_t nu2 = val[2]*val[2];
1373  Double_t Gamma_t = 2.0*val[1]*val[1]/((wL2+nu2)*(wL2+nu2))*
1374  ((wL2+nu2)*val[2]*t
1375  + (wL2-nu2)*(1.0 - TMath::Exp(-val[2]*t)*TMath::Cos(wL*t))
1376  - 2.0*val[2]*wL*TMath::Exp(-val[2]*t)*TMath::Sin(wL*t));
1377  result = TMath::Exp(-Gamma_t);
1378  } else { // from Voltera
1379  result = GetDynKTLFValue(tt);
1380  }
1381 
1382  return result;
1383 
1384 }
1385 
1386 //--------------------------------------------------------------------------
1401 Double_t PTheory::StaticLorentzKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1402 {
1403  // expected parameters: lambda [tshift]
1404 
1405  Double_t val[2];
1406 
1407  assert(fParamNo.size() <= 2);
1408 
1409  // check if FUNCTIONS are used
1410  for (UInt_t i=0; i<fParamNo.size(); i++) {
1411  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1412  val[i] = paramValues[fParamNo[i]];
1413  } else { // function
1414  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1415  }
1416  }
1417 
1418  Double_t a_t;
1419  if (fParamNo.size() == 1) // no tshift
1420  a_t = t*val[0];
1421  else // tshift present
1422  a_t = (t-val[1])*val[0];
1423 
1424  return 0.333333333333333 * (1.0 + 2.0*(1.0 - a_t)*TMath::Exp(-a_t));
1425 }
1426 
1427 //--------------------------------------------------------------------------
1444 Double_t PTheory::StaticLorentzKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1445 {
1446  // expected parameters: frequency damping [tshift]
1447 
1448  Double_t val[3];
1449  Double_t result;
1450 
1451  assert(fParamNo.size() <= 3);
1452 
1453  // check if FUNCTIONS are used
1454  for (UInt_t i=0; i<fParamNo.size(); i++) {
1455  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1456  val[i] = paramValues[fParamNo[i]];
1457  } else { // function
1458  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1459  }
1460  }
1461 
1462  // check if all parameters == 0
1463  if ((val[0] == 0.0) && (val[1] == 0.0))
1464  return 1.0;
1465 
1466  // check if the parameter values have changed, and if yes recalculate the non-analytic integral
1467  // check only the first two parameters since the tshift is irrelevant for the LF-integral calculation!!
1468  Bool_t newParam = false;
1469  for (UInt_t i=0; i<2; i++) {
1470  if (val[i] != fPrevParam[i]) {
1471  newParam = true;
1472  break;
1473  }
1474  }
1475 
1476  if (newParam) { // new parameters found
1477  for (UInt_t i=0; i<2; i++)
1478  fPrevParam[i] = val[i];
1479  CalculateLorentzLFIntegral(val);
1480  }
1481 
1482  Double_t tt;
1483  if (fParamNo.size() == 2) // no tshift
1484  tt = t;
1485  else // tshift present
1486  tt = t-val[2];
1487 
1488  if (tt < 0.0) // for times < 0 return a function value of 1.0
1489  return 1.0;
1490 
1491  if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use the ZF formula
1492  Double_t at = tt*val[1];
1493  result = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
1494  } else if (val[1]/val[0] > 159.1549) { // check if a/w0 > 1000.0, in which case the ZF formula is used
1495  Double_t at = tt*val[1];
1496  result = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
1497  } else {
1498  Double_t a = val[1];
1499  Double_t at = a*tt;
1500  Double_t w0 = 2.0*TMath::Pi()*val[0];
1501  Double_t a_w0 = a/w0;
1502  Double_t w0t = w0*tt;
1503 
1504  Double_t j1, j0;
1505  if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
1506  j0 = 1.0;
1507  j1 = 0.0;
1508  } else {
1509  j0 = sin(w0t)/w0t;
1510  j1 = (sin(w0t)-w0t*cos(w0t))/(w0t*w0t);
1511  }
1512 
1513  result = 1.0 - a_w0*j1*exp(-at) - a_w0*a_w0*(j0*exp(-at) - 1.0) - GetLFIntegralValue(tt);
1514  }
1515 
1516  return result;
1517 
1518 }
1519 
1520 //--------------------------------------------------------------------------
1541 Double_t PTheory::DynamicLorentzKTLF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1542 {
1543  // expected parameters: frequency damping hopping [tshift]
1544 
1545  Double_t val[4];
1546  Double_t result = 0.0;
1547 
1548  assert(fParamNo.size() <= 4);
1549 
1550  // check if FUNCTIONS are used
1551  for (UInt_t i=0; i<fParamNo.size(); i++) {
1552  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1553  val[i] = paramValues[fParamNo[i]];
1554  } else { // function
1555  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1556  }
1557  }
1558 
1559  // check if all parameters == 0
1560  if ((val[0] == 0.0) && (val[1] == 0.0) && (val[2] == 0.0))
1561  return 1.0;
1562 
1563  // make sure that damping and hopping are positive definite
1564  if (val[1] < 0.0)
1565  val[1] = -val[1];
1566  if (val[2] < 0.0)
1567  val[2] = -val[2];
1568 
1569 
1570  Double_t tt;
1571  if (fParamNo.size() == 3) // no tshift
1572  tt = t;
1573  else // tshift present
1574  tt = t-val[3];
1575 
1576  if (tt < 0.0) // for times < 0 return a function value of 1.0
1577  return 1.0;
1578 
1579  // check if hopping > 5 * damping, of Larmor angular frequency is > 30 * damping (BMW limit)
1580  Double_t w0 = 2.0*TMath::Pi()*val[0];
1581  Double_t a = val[1];
1582  Double_t nu = val[2];
1583  if ((nu > 5.0 * a) || (w0 >= 30.0 * a)) {
1584  // 'c' and 'd' are parameters BMW obtained by fitting large parameter space LF-curves to the model below
1585  const Double_t c[7] = {1.15331, 1.64826, -0.71763, 3.0, 0.386683, -5.01876, 2.41854};
1586  const Double_t d[4] = {2.44056, 2.92063, 1.69581, 0.667277};
1587  Double_t w0N[4];
1588  Double_t nuN[4];
1589  w0N[0] = w0;
1590  nuN[0] = nu;
1591  for (UInt_t i=1; i<4; i++) {
1592  w0N[i] = w0 * w0N[i-1];
1593  nuN[i] = nu * nuN[i-1];
1594  }
1595  Double_t denom = w0N[3]+d[0]*w0N[2]*nuN[0]+d[1]*w0N[1]*nuN[1]+d[2]*w0N[0]*nuN[2]+d[3]*nuN[3];
1596  Double_t b1 = (c[0]*w0N[2]+c[1]*w0N[1]*nuN[0]+c[2]*w0N[0]*nuN[1])/denom;
1597  Double_t b2 = (c[3]*w0N[2]+c[4]*w0N[1]*nuN[0]+c[5]*w0N[0]*nuN[1]+c[6]*nuN[2])/denom;
1598 
1599  Double_t w0t = w0*tt;
1600  Double_t j1, j0;
1601  if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
1602  j0 = 1.0;
1603  j1 = 0.0;
1604  } else {
1605  j0 = sin(w0t)/w0t;
1606  j1 = (sin(w0t)-w0t*cos(w0t))/(w0t*w0t);
1607  }
1608 
1609  Double_t Gamma_t = -4.0/3.0*a*(b1*(1.0-j0*TMath::Exp(-nu*tt))+b2*j1*TMath::Exp(-nu*tt)+(1.0-b2*w0/3.0-b1*nu)*tt);
1610 
1611  return TMath::Exp(Gamma_t);
1612  }
1613 
1614  // check if the parameter values have changed, and if yes recalculate the non-analytic integral
1615  // check only the first three parameters since the tshift is irrelevant for the LF-integral calculation!!
1616  Bool_t newParam = false;
1617  for (UInt_t i=0; i<3; i++) {
1618  if (val[i] != fPrevParam[i]) {
1619  newParam = true;
1620  break;
1621  }
1622  }
1623 
1624  if (newParam) { // new parameters found
1625  for (UInt_t i=0; i<3; i++)
1626  fPrevParam[i] = val[i];
1627  CalculateDynKTLF(val, 1); // 1 means Lorentz
1628  }
1629 
1630  result = GetDynKTLFValue(tt);
1631 
1632  return result;
1633 
1634 }
1635 
1636 //--------------------------------------------------------------------------
1651 Double_t PTheory::CombiLGKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1652 {
1653  // expected parameters: lambdaL lambdaG [tshift]
1654 
1655  Double_t val[3];
1656 
1657  assert(fParamNo.size() <= 3);
1658 
1659  // check if FUNCTIONS are used
1660  for (UInt_t i=0; i<fParamNo.size(); i++) {
1661  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1662  val[i] = paramValues[fParamNo[i]];
1663  } else { // function
1664  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1665  }
1666  }
1667 
1668  Double_t tt;
1669  if (fParamNo.size() == 2) // no tshift
1670  tt = t;
1671  else // tshift present
1672  tt = t-val[2];
1673 
1674  Double_t lambdaL_t = tt*val[0];
1675  Double_t lambdaG_t_2 = tt*tt*val[1]*val[1];
1676 
1677  return 0.333333333333333 *
1678  (1.0 + 2.0*(1.0-lambdaL_t-lambdaG_t_2)*TMath::Exp(-(lambdaL_t+0.5*lambdaG_t_2)));
1679 }
1680 
1681 //--------------------------------------------------------------------------
1697 Double_t PTheory::StrKT(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1698 {
1699  // expected parameters: sigma beta [tshift]
1700 
1701  Double_t val[3];
1702 
1703  assert(fParamNo.size() <= 3);
1704 
1705  // check if FUNCTIONS are used
1706  for (UInt_t i=0; i<fParamNo.size(); i++) {
1707  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1708  val[i] = paramValues[fParamNo[i]];
1709  } else { // function
1710  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1711  }
1712  }
1713 
1714  // check for beta too small (beta < 0.1) in which case numerical problems could arise and the function is anyhow
1715  // almost identical to a constant of 1/3.
1716  if (val[1] < 0.1)
1717  return 0.333333333333333;
1718 
1719  Double_t tt;
1720  if (fParamNo.size() == 2) // no tshift
1721  tt = t;
1722  else // tshift present
1723  tt = t-val[2];
1724 
1725  Double_t sigma_t = TMath::Power(tt*val[0],val[1]);
1726 
1727  return 0.333333333333333 *
1728  (1.0 + 2.0*(1.0-sigma_t)*TMath::Exp(-sigma_t/val[1]));
1729 }
1730 
1731 //--------------------------------------------------------------------------
1746 Double_t PTheory::SpinGlass(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1747 {
1748  // expected parameters: lambda gamma q [tshift]
1749 
1750  if (paramValues[fParamNo[0]] == 0.0)
1751  return 1.0;
1752 
1753  Double_t val[4];
1754 
1755  assert(fParamNo.size() <= 4);
1756 
1757  // check if FUNCTIONS are used
1758  for (UInt_t i=0; i<fParamNo.size(); i++) {
1759  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1760  val[i] = paramValues[fParamNo[i]];
1761  } else { // function
1762  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1763  }
1764  }
1765 
1766  Double_t tt;
1767  if (fParamNo.size() == 3) // no tshift
1768  tt = t;
1769  else // tshift present
1770  tt = t-val[3];
1771 
1772  Double_t lambda_2 = val[0]*val[0];
1773  Double_t lambda_t_2_q = tt*tt*lambda_2*val[2];
1774  Double_t rate_2 = 4.0*lambda_2*(1.0-val[2])*tt/val[1];
1775 
1776  Double_t rateL = TMath::Sqrt(fabs(rate_2));
1777  Double_t rateT = TMath::Sqrt(fabs(rate_2)+lambda_t_2_q);
1778 
1779  return 0.333333333333333*(TMath::Exp(-rateL) + 2.0*(1.0-lambda_t_2_q/rateT)*TMath::Exp(-rateT));
1780 }
1781 
1782 //--------------------------------------------------------------------------
1797 Double_t PTheory::RandomAnisotropicHyperfine(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1798 {
1799  // expected parameters: nu lambda [tshift]
1800 
1801  Double_t val[3];
1802 
1803  assert(fParamNo.size() <= 3);
1804 
1805  // check if FUNCTIONS are used
1806  for (UInt_t i=0; i<fParamNo.size(); i++) {
1807  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1808  val[i] = paramValues[fParamNo[i]];
1809  } else { // function
1810  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1811  }
1812  }
1813 
1814  Double_t tt;
1815  if (fParamNo.size() == 2) // no tshift
1816  tt = t;
1817  else // tshift present
1818  tt = t-val[2];
1819 
1820  Double_t nu_t = tt*val[0];
1821  Double_t lambda_t = tt*val[1];
1822 
1823  return 0.166666666666667*(1.0-0.5*nu_t)*TMath::Exp(-0.5*nu_t) +
1824  0.333333333333333*(1.0-0.25*nu_t)*TMath::Exp(-0.25*(nu_t+2.44949*lambda_t));
1825 }
1826 
1827 //--------------------------------------------------------------------------
1842 Double_t PTheory::Abragam(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1843 {
1844  // expected parameters: sigma gamma [tshift]
1845 
1846  Double_t val[3];
1847 
1848  assert(fParamNo.size() <= 3);
1849 
1850  // check if FUNCTIONS are used
1851  for (UInt_t i=0; i<fParamNo.size(); i++) {
1852  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1853  val[i] = paramValues[fParamNo[i]];
1854  } else { // function
1855  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1856  }
1857  }
1858 
1859  Double_t tt;
1860  if (fParamNo.size() == 2) // no tshift
1861  tt = t;
1862  else // tshift present
1863  tt = t-val[2];
1864 
1865  Double_t gamma_t = tt*val[1];
1866 
1867  return TMath::Exp(-TMath::Power(val[0]/val[1],2.0)*
1868  (TMath::Exp(-gamma_t)-1.0+gamma_t));
1869 }
1870 
1871 //--------------------------------------------------------------------------
1886 Double_t PTheory::TFCos(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1887 {
1888  // expected parameters: phase frequency [tshift]
1889 
1890  Double_t val[3];
1891 
1892  assert(fParamNo.size() <= 3);
1893 
1894  // check if FUNCTIONS are used
1895  for (UInt_t i=0; i<fParamNo.size(); i++) {
1896  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1897  val[i] = paramValues[fParamNo[i]];
1898  } else { // function
1899  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1900  }
1901  }
1902 
1903  Double_t tt;
1904  if (fParamNo.size() == 2) // no tshift
1905  tt = t;
1906  else // tshift present
1907  tt = t-val[2];
1908 
1909  return TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
1910 }
1911 
1912 //--------------------------------------------------------------------------
1927 Double_t PTheory::InternalField(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1928 {
1929  // expected parameters: fraction phase frequency rateT rateL [tshift]
1930 
1931  Double_t val[6];
1932 
1933  assert(fParamNo.size() <= 6);
1934 
1935  // check if FUNCTIONS are used
1936  for (UInt_t i=0; i<fParamNo.size(); i++) {
1937  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1938  val[i] = paramValues[fParamNo[i]];
1939  } else { // function
1940  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1941  }
1942  }
1943 
1944  Double_t tt;
1945  if (fParamNo.size() == 5) // no tshift
1946  tt = t;
1947  else // tshift present
1948  tt = t-val[5];
1949 
1950  return val[0]*TMath::Cos(DEG_TO_RAD*val[1]+TWO_PI*val[2]*tt)*TMath::Exp(-val[3]*tt) +
1951  (1-val[0])*TMath::Exp(-val[4]*tt);
1952 }
1953 
1954 //--------------------------------------------------------------------------
1969 Double_t PTheory::InternalFieldGK(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
1970 {
1971  // expected parameters: [0]:fraction [1]:frequency [2]:sigma [3]:lambda [4]:beta [[5]:tshift]
1972 
1973  Double_t val[6];
1974 
1975  assert(fParamNo.size() <= 6);
1976 
1977  // check if FUNCTIONS are used
1978  for (UInt_t i=0; i<fParamNo.size(); i++) {
1979  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
1980  val[i] = paramValues[fParamNo[i]];
1981  } else { // function
1982  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
1983  }
1984  }
1985 
1986  Double_t tt;
1987  if (fParamNo.size() == 5) // no tshift
1988  tt = t;
1989  else // tshift present
1990  tt = t-val[5];
1991 
1992  Double_t result = 0.0;
1993  Double_t w_t = TWO_PI*val[1]*tt;
1994  Double_t rateLF = TMath::Power(val[3]*tt, val[4]);
1995  Double_t rate2 = val[2]*val[2]*tt*tt; // (sigma t)^2
1996 
1997  if (val[1] < 0.01) { // internal field frequency is approaching zero
1998  result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(1.0-rate2)*TMath::Exp(-0.5*rate2);
1999  } else {
2000  result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(TMath::Cos(w_t)-val[2]*val[2]*tt/(TWO_PI*val[1])*TMath::Sin(w_t))*TMath::Exp(-0.5*rate2);
2001  }
2002 
2003  return result;
2004 }
2005 
2006 //--------------------------------------------------------------------------
2021 Double_t PTheory::InternalFieldLL(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2022 {
2023  // expected parameters: [0]:fraction [1]:frequency [2]:a [3]:lambda [4]:beta [[5]:tshift]
2024 
2025  Double_t val[6];
2026 
2027  assert(fParamNo.size() <= 6);
2028 
2029  // check if FUNCTIONS are used
2030  for (UInt_t i=0; i<fParamNo.size(); i++) {
2031  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2032  val[i] = paramValues[fParamNo[i]];
2033  } else { // function
2034  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2035  }
2036  }
2037 
2038  Double_t tt;
2039  if (fParamNo.size() == 5) // no tshift
2040  tt = t;
2041  else // tshift present
2042  tt = t-val[5];
2043 
2044  Double_t result = 0.0;
2045  Double_t w_t = TWO_PI*val[1]*tt;
2046  Double_t rateLF = TMath::Power(val[3]*tt, val[4]);
2047  Double_t a_t = val[2]*tt; // a t
2048 
2049  if (val[1] < 0.01) { // internal field frequency is approaching zero
2050  result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(1.0-a_t)*TMath::Exp(-a_t);
2051  } else {
2052  result = (1.0-val[0])*TMath::Exp(-rateLF) + val[0]*(TMath::Cos(w_t)-val[3]/(TWO_PI*val[1])*TMath::Sin(w_t))*TMath::Exp(-a_t);
2053  }
2054 
2055  return result;
2056 }
2057 
2058 //--------------------------------------------------------------------------
2073 Double_t PTheory::Bessel(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2074 {
2075  // expected parameters: phase frequency [tshift]
2076 
2077  Double_t val[3];
2078 
2079  assert(fParamNo.size() <= 3);
2080 
2081  // check if FUNCTIONS are used
2082  for (UInt_t i=0; i<fParamNo.size(); i++) {
2083  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2084  val[i] = paramValues[fParamNo[i]];
2085  } else { // function
2086  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2087  }
2088  }
2089 
2090  Double_t tt;
2091  if (fParamNo.size() == 2) // no tshift
2092  tt = t;
2093  else // tshift present
2094  tt = t-val[2];
2095 
2096  return TMath::BesselJ0(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
2097 }
2098 
2099 //--------------------------------------------------------------------------
2114 Double_t PTheory::InternalBessel(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2115 {
2116  // expected parameters: fraction phase frequency rateT rateL [tshift]
2117 
2118  Double_t val[6];
2119 
2120  assert(fParamNo.size() <= 6);
2121 
2122  // check if FUNCTIONS are used
2123  for (UInt_t i=0; i<fParamNo.size(); i++) {
2124  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2125  val[i] = paramValues[fParamNo[i]];
2126  } else { // function
2127  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2128  }
2129  }
2130 
2131  Double_t tt;
2132  if (fParamNo.size() == 5) // no tshift
2133  tt = t;
2134  else // tshift present
2135  tt = t-val[5];
2136 
2137  return val[0]* TMath::BesselJ0(DEG_TO_RAD*val[1]+TWO_PI*val[2]*tt)*
2138  TMath::Exp(-val[3]*tt) +
2139  (1.0-val[0])*TMath::Exp(-val[4]*tt);
2140 }
2141 
2142 //--------------------------------------------------------------------------
2163 Double_t PTheory::SkewedGauss(Double_t t, const PDoubleVector &paramValues,
2164  const PDoubleVector &funcValues) const {
2165  // Expected parameters: phase, frequency, sigma-, sigma+, [tshift].
2166  // To be stored in the array "val" as:
2167  // val[0] = phase
2168  // val[1] = frequency
2169  // val[2] = sigma-
2170  // val[3] = sigma+
2171  // val[4] = tshift [optional]
2172  Double_t val[5];
2173 
2174  // Check that we have the correct number of fit parameters.
2175  assert(fParamNo.size() <= 5);
2176 
2177  // Check if FUNCTIONS are used.
2178  for (UInt_t i = 0; i < fParamNo.size(); i++) {
2179  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2180  val[i] = paramValues[fParamNo[i]];
2181  } else { // function
2182  val[i] = funcValues[fParamNo[i] - MSR_PARAM_FUN_OFFSET];
2183  }
2184  }
2185 
2186  // Apply the tshift (if required).
2187  Double_t tt = t;
2188  if (fParamNo.size() == 5) {
2189  tt = t - val[4];
2190  }
2191 
2192  // Evaluate the skewed Gaussian!
2193 
2194  // First, calculate some "helper" terms.
2195  Double_t sigma_p = std::abs(val[3]);
2196  Double_t sigma_m = std::abs(val[2]);
2197  Double_t arg_p = sigma_p * tt;
2198  Double_t arg_m = sigma_m * tt;
2199  Double_t z_p = arg_p / SQRT_TWO; // sigma+
2200  Double_t z_m = arg_m / SQRT_TWO; // sigma-
2201  Double_t g_p = TMath::Exp(-0.5 * arg_p * arg_p); // gauss sigma+
2202  Double_t g_m = TMath::Exp(-0.5 * arg_m * arg_m); // gauss sigma-
2203  Double_t w_p = sigma_p / (sigma_p + sigma_m);
2204  Double_t w_m = 1.0 - w_p;
2205  Double_t phase = DEG_TO_RAD * val[0];
2206  Double_t freq = TWO_PI * val[1];
2207 
2208  // Evalute the EVEN frequency component of the skewed Gaussian.
2209  Double_t skg_cos = TMath::Cos(phase + freq * tt) * (w_m * g_m + w_p * g_p);
2210 
2211  // Evalute the ODD frequency component of the skewed Gaussian.
2212  constexpr Double_t z_max = 26.7776;
2213  // Note: the check against z_max is needed to prevent floating-point overflow
2214  // in the return value of ROOT::Math::conf_hyperg(1/2, 3/2, z * z)
2215  // (i.e., confluent hypergeometric function of the first kind, 1F1).
2216  // In the case that z > z_max, return zero (otherwise there is some
2217  // numeric discontinuity at later times).
2218  Double_t skg_sin =
2219  TMath::Sin(phase + freq * tt) *
2220  ((z_m > z_max) or (z_p > z_max)
2221  ? 0.0
2222  : (w_m * g_m * 2.0 * z_m / SQRT_PI) *
2223  ROOT::Math::conf_hyperg(0.5, 1.5, z_m * z_m) -
2224  (w_p * g_p * 2.0 * z_p / SQRT_PI) *
2225  ROOT::Math::conf_hyperg(0.5, 1.5, z_p * z_p));
2226 
2227  // Return the skewed Gaussian: skg = skg_cos + skg_sin.
2228  // Also check that skg_sin is finite!
2229  return skg_cos + (std::isfinite(skg_sin) ? skg_sin : 0.0);
2230 }
2231 
2232 //--------------------------------------------------------------------------
2252 Double_t PTheory::StaticNKZF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2253 {
2254  // expected paramters: damping_D0 [0], R_b tshift [1]
2255 
2256  Double_t val[3];
2257  Double_t result = 1.0;
2258 
2259  assert(fParamNo.size() <= 3);
2260 
2261  if (t < 0.0)
2262  return result;
2263 
2264  // check if FUNCTIONS are used
2265  for (UInt_t i=0; i<fParamNo.size(); i++) {
2266  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2267  val[i] = paramValues[fParamNo[i]];
2268  } else { // function
2269  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2270  }
2271  }
2272 
2273  Double_t tt;
2274  if (fParamNo.size() == 2) // no tshift
2275  tt = t;
2276  else // tshift present
2277  tt = t-val[2];
2278 
2279  Double_t D2_t2 = val[0]*val[0]*tt*tt;
2280  Double_t denom = 1.0+val[1]*val[1]*D2_t2;
2281 
2282  result = 0.333333333333333 + 0.666666666666666667 * TMath::Power(1.0/denom, 1.5) * (1.0 - (D2_t2/denom)) * exp(-0.5*D2_t2/denom);
2283 
2284  return result;
2285 }
2286 
2287 //--------------------------------------------------------------------------
2307 Double_t PTheory::StaticNKTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2308 {
2309  // expected paramters: phase [0], frequency [1], damping_D0 [2], R_b [3], [tshift [4]]
2310 
2311  Double_t val[5];
2312  Double_t result = 1.0;
2313 
2314  assert(fParamNo.size() <= 5);
2315 
2316  if (t < 0.0)
2317  return result;
2318 
2319  // check if FUNCTIONS are used
2320  for (UInt_t i=0; i<fParamNo.size(); i++) {
2321  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2322  val[i] = paramValues[fParamNo[i]];
2323  } else { // function
2324  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2325  }
2326  }
2327 
2328  Double_t tt;
2329  if (fParamNo.size() == 4) // no tshift
2330  tt = t;
2331  else // tshift present
2332  tt = t-val[4];
2333 
2334  Double_t D2_t2 = val[2]*val[2]*tt*tt;
2335  Double_t denom = 1.0+val[3]*val[3]*D2_t2;
2336 
2337  result = sqrt(1.0/denom)*exp(-0.5*D2_t2/denom)*TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
2338 
2339  return result;
2340 }
2341 
2342 //--------------------------------------------------------------------------
2363 Double_t PTheory::DynamicNKZF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2364 {
2365  // expected paramters: damping_D0 [0], R_b [1], nu_c [2], [tshift [3]]
2366 
2367  Double_t val[4];
2368  Double_t result = 1.0;
2369 
2370  assert(fParamNo.size() <= 4);
2371 
2372  if (t < 0.0)
2373  return result;
2374 
2375  // check if FUNCTIONS are used
2376  for (UInt_t i=0; i<fParamNo.size(); i++) {
2377  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2378  val[i] = paramValues[fParamNo[i]];
2379  } else { // function
2380  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2381  }
2382  }
2383 
2384  Double_t tt;
2385  if (fParamNo.size() == 3) // no tshift
2386  tt = t;
2387  else // tshift present
2388  tt = t-val[3];
2389 
2390  Double_t theta;
2391  if (val[2] < 1.0e-6) { // nu_c -> 0
2392  theta = 0.5*tt*tt;
2393  } else {
2394  theta = (exp(-val[2]*tt) - 1.0 + val[2]*tt)/(val[2]*val[2]);
2395  }
2396  Double_t denom = 1.0 + 4.0*val[0]*val[0]*val[1]*val[1]*theta;
2397 
2398  result = sqrt(1.0/denom)*exp(-2.0*val[0]*val[0]*theta/denom);
2399 
2400  return result;
2401 }
2402 
2403 //--------------------------------------------------------------------------
2424 Double_t PTheory::DynamicNKTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2425 {
2426  // expected paramters: phase [0], frequency [1], damping_D0 [2], R_b [3], nu_c [4], [tshift [5]]
2427 
2428  Double_t val[6];
2429  Double_t result = 1.0;
2430 
2431  assert(fParamNo.size() <= 6);
2432 
2433  if (t < 0.0)
2434  return result;
2435 
2436  // check if FUNCTIONS are used
2437  for (UInt_t i=0; i<fParamNo.size(); i++) {
2438  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2439  val[i] = paramValues[fParamNo[i]];
2440  } else { // function
2441  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2442  }
2443  }
2444 
2445  Double_t tt;
2446  if (fParamNo.size() == 5) // no tshift
2447  tt = t;
2448  else // tshift present
2449  tt = t-val[5];
2450 
2451  Double_t theta;
2452  if (val[4] < 1.0e-6) { // nu_c -> 0
2453  theta = 0.5*tt*tt;
2454  } else {
2455  theta = (exp(-val[4]*tt) - 1.0 + val[4]*tt)/(val[4]*val[4]);
2456  }
2457  Double_t denom = 1.0 + 2.0*val[2]*val[2]*val[3]*val[3]*theta;
2458 
2459  result = sqrt(1.0/denom)*exp(-val[2]*val[2]*theta/denom)*TMath::Cos(DEG_TO_RAD*val[0]+TWO_PI*val[1]*tt);
2460 
2461  return result;
2462 }
2463 
2464 //--------------------------------------------------------------------------
2478 Double_t PTheory::Polynom(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2479 {
2480  // expected parameters: tshift p0 p1 p2 ...
2481 
2482  Double_t result = 0.0;
2483  Double_t tshift = 0.0;
2484  Double_t val;
2485  Double_t expo = 0.0;
2486 
2487  // check if FUNCTIONS are used
2488  for (UInt_t i=0; i<fParamNo.size(); i++) {
2489  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2490  val = paramValues[fParamNo[i]];
2491  } else { // function
2492  val = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2493  }
2494  if (i==0) { // tshift
2495  tshift = val;
2496  continue;
2497  }
2498  result += val*pow(t-tshift, expo);
2499  expo++;
2500  }
2501 
2502  return result;
2503 }
2504 
2505 //--------------------------------------------------------------------------
2515 Double_t PTheory::UserFcn(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2516 {
2517  // check if FUNCTIONS are used
2518  for (UInt_t i=0; i<fUserParam.size(); i++) {
2519  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2520  fUserParam[i] = paramValues[fParamNo[i]];
2521  } else { // function
2522  fUserParam[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2523  }
2524  }
2525 
2526  return (*fUserFcn)(t, fUserParam);
2527 }
2528 
2529 //--------------------------------------------------------------------------
2542 void PTheory::CalculateGaussLFIntegral(const Double_t *val) const
2543 {
2544  // val[0] = nu (field), val[1] = Delta
2545 
2546  if (val[0] == 0.0) { // field == 0.0, hence nothing to be done
2547  return;
2548  } else if (val[1]/val[0] > 79.5775) { // check if a/w0 > 500.0, in which case the ZF formula is used and here nothing has to be done
2549  return;
2550  }
2551 
2552 
2553  Double_t dt=0.001; // all times in usec
2554  Double_t t, ft;
2555  Double_t w0 = TMath::TwoPi()*val[0];
2556  Double_t Delta = val[1];
2557  Double_t preFactor = 2.0*TMath::Power(Delta, 4.0) / TMath::Power(w0, 3.0);
2558 
2559  // check if the time resolution needs to be increased
2560  const Int_t samplingPerPeriod = 20;
2561  const Int_t samplingOnExp = 3000;
2562  if ((Delta <= w0) && (1.0/val[0] < 20.0)) { // makes sure that the frequency sampling is fine enough
2563  if (1.0/val[0]/samplingPerPeriod < 0.001) {
2564  dt = 1.0/val[0]/samplingPerPeriod;
2565  }
2566  } else if ((Delta > w0) && (Delta <= 10.0)) {
2567  if (Delta/w0 > 10.0) {
2568  dt = 0.00005;
2569  }
2570  } else if ((Delta > w0) && (Delta > 10.0)) { // makes sure there is a fine enough sampling for large Delta's
2571  if (1.0/Delta/samplingOnExp < 0.001) {
2572  dt = 1.0/Delta/samplingOnExp;
2573  }
2574  }
2575 
2576  // keep sampling time
2577  fSamplingTime = dt;
2578 
2579  // clear previously allocated vector
2580  fLFIntegral.clear();
2581 
2582  // calculate integral
2583  t = 0.0;
2584  fLFIntegral.push_back(0.0); // start value of the integral
2585 
2586  ft = 0.0;
2587  Double_t step = 0.0, lastft = 1.0, diff = 0.0;
2588  do {
2589  t += dt;
2590  step = 0.5*dt*preFactor*(exp(-0.5*pow(Delta * (t-dt), 2.0))*sin(w0*(t-dt))+
2591  exp(-0.5*pow(Delta * t, 2.0))*sin(w0*t));
2592  ft += step;
2593  diff = fabs(fabs(lastft)-fabs(ft));
2594  lastft = ft;
2595  fLFIntegral.push_back(ft);
2596  } while ((t <= 20.0) && (diff > 1.0e-10));
2597 }
2598 
2599 //--------------------------------------------------------------------------
2612 void PTheory::CalculateLorentzLFIntegral(const Double_t *val) const
2613 {
2614  // val[0] = nu, val[1] = a
2615 
2616  // a few checks if the integral actually needs to be calculated
2617  if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use the ZF formula and here nothing has to be done
2618  return;
2619  } else if (val[1]/val[0] > 159.1549) { // check if a/w0 > 1000.0, in which case the ZF formula is used and here nothing has to be done
2620  return;
2621  }
2622 
2623  Double_t dt=0.001; // all times in usec
2624  Double_t t, ft;
2625  Double_t w0 = TMath::TwoPi()*val[0];
2626  Double_t a = val[1];
2627  Double_t preFactor = a*(1+pow(a/w0,2.0));
2628 
2629  // check if the time resolution needs to be increased
2630  const Int_t samplingPerPeriod = 20;
2631  const Int_t samplingOnExp = 3000;
2632  if ((a <= w0) && (1.0/val[0] < 20.0)) { // makes sure that the frequency sampling is fine enough
2633  if (1.0/val[0]/samplingPerPeriod < 0.001) {
2634  dt = 1.0/val[0]/samplingPerPeriod;
2635  }
2636  } else if ((a > w0) && (a <= 10.0)) {
2637  if (a/w0 > 10.0) {
2638  dt = 0.00005;
2639  }
2640  } else if ((a > w0) && (a > 10.0)) { // makes sure there is a fine enough sampling for large a's
2641  if (1.0/a/samplingOnExp < 0.001) {
2642  dt = 1.0/a/samplingOnExp;
2643  }
2644  }
2645 
2646  // keep sampling time
2647  fSamplingTime = dt;
2648 
2649  // clear previously allocated vector
2650  fLFIntegral.clear();
2651 
2652  // calculate integral
2653  t = 0.0;
2654  fLFIntegral.push_back(0.0); // start value of the integral
2655 
2656  ft = 0.0;
2657  // calculate first integral bin value (needed bcause of sin(x)/x x->0)
2658  t += dt;
2659  ft += 0.5*dt*preFactor*(1.0+sin(w0*t)/(w0*t)*exp(-a*t));
2660  fLFIntegral.push_back(ft);
2661  // calculate all the other integral bin values
2662  Double_t step = 0.0, lastft = 1.0, diff = 0.0;
2663  do {
2664  t += dt;
2665  step = 0.5*dt*preFactor*(sin(w0*(t-dt))/(w0*(t-dt))*exp(-a*(t-dt))+sin(w0*t)/(w0*t)*exp(-a*t));
2666  ft += step;
2667  diff = fabs(fabs(lastft)-fabs(ft));
2668  lastft = ft;
2669  fLFIntegral.push_back(ft);
2670  } while ((t <= 20.0) && (diff > 1.0e-10));
2671 }
2672 
2673 
2674 //--------------------------------------------------------------------------
2682 Double_t PTheory::GetLFIntegralValue(const Double_t t) const
2683 {
2684  if (t < 0.0)
2685  return 0.0;
2686 
2687  UInt_t idx = static_cast<UInt_t>(t/fSamplingTime);
2688 
2689  if (idx + 2 > fLFIntegral.size())
2690  return fLFIntegral.back();
2691 
2692  // linearly interpolate between the two relevant function bins
2693  Double_t df = (fLFIntegral[idx+1]-fLFIntegral[idx])*(t/fSamplingTime-static_cast<Double_t>(idx));
2694 
2695  return fLFIntegral[idx]+df;
2696 }
2697 
2698 //--------------------------------------------------------------------------
2709 void PTheory::CalculateDynKTLF(const Double_t *val, Int_t tag) const
2710 {
2711  // val: 0=nu0, 1=Delta (Gauss) / a (Lorentz), 2=nu
2712  const Double_t Tmax = 20.0; // 20 usec
2713  UInt_t N = static_cast<UInt_t>(16.0*Tmax*val[0]);
2714 
2715  // check if rate (Delta or a) is very high
2716  if (fabs(val[1]) > 0.1) {
2717  Double_t tmin = 20.0;
2718  switch (tag) {
2719  case 0: // Gauss
2720  tmin = fabs(sqrt(3.0)/val[1]);
2721  break;
2722  case 1: // Lorentz
2723  tmin = fabs(2.0/val[1]);
2724  break;
2725  default:
2726  break;
2727  }
2728  UInt_t Nrate = static_cast<UInt_t>(25.0 * Tmax / tmin);
2729  if (Nrate > N) {
2730  N = Nrate;
2731  }
2732  }
2733 
2734  if (N < 300) // if too few points, i.e. nu0 very small, take 300 points
2735  N = 300;
2736 
2737  if (N>1e6) // make sure that N is not too large to prevent memory overflow
2738  N = 1e6;
2739 
2740  // allocate memory for dyn KT LF function vector
2741  fDynLFFuncValue.clear(); // get rid of a possible previous vector
2742  fDynLFFuncValue.resize(N);
2743 
2744  // calculate the non-analytic integral of the static KT LF function
2745  switch (tag) {
2746  case 0: // Gauss
2747  CalculateGaussLFIntegral(val);
2748  break;
2749  case 1: // Lorentz
2750  CalculateLorentzLFIntegral(val);
2751  break;
2752  default:
2753  std::cerr << std::endl << ">> PTheory::CalculateDynKTLF: **FATAL ERROR** You should never have reached this point." << std::endl;
2754  assert(false);
2755  break;
2756  }
2757 
2758  // calculate the P^(0)(t) exp(-nu t) vector
2759  PDoubleVector p0exp(N);
2760  Double_t t = 0.0;
2761  Double_t dt = Tmax/N;
2762  fDynLFdt = dt; // keep it since it is needed in GetDynKTLFValue()
2763  for (UInt_t i=0; i<N; i++) {
2764  switch (tag) {
2765  case 0: // Gauss
2766  if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use zero field formula
2767  Double_t sigma_t_2 = t*t*val[1]*val[1];
2768  p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
2769  } else if (val[1]/val[0] > 79.5775) { // check if Delta/w0 > 500.0, in which case the ZF formula is used
2770  Double_t sigma_t_2 = t*t*val[1]*val[1];
2771  p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - sigma_t_2)*TMath::Exp(-0.5*sigma_t_2));
2772  } else {
2773  Double_t delta = val[1];
2774  Double_t w0 = TWO_PI*val[0];
2775 
2776  p0exp[i] = 1.0 - 2.0*TMath::Power(delta/w0,2.0)*(1.0 -
2777  TMath::Exp(-0.5*TMath::Power(delta*t, 2.0))*TMath::Cos(w0*t)) +
2778  GetLFIntegralValue(t);
2779  }
2780  break;
2781  case 1: // Lorentz
2782  if (val[0] < 0.02) { // if smaller 20kHz ~ 0.27G use zero field formula
2783  Double_t at = t*val[1];
2784  p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
2785  } else if (val[1]/val[0] > 159.1549) { // check if a/w0 > 1000.0, in which case the ZF formula is used
2786  Double_t at = t*val[1];
2787  p0exp[i] = 0.333333333333333 * (1.0 + 2.0*(1.0 - at)*TMath::Exp(-at));
2788  } else {
2789  Double_t a = val[1];
2790  Double_t at = a*t;
2791  Double_t w0 = TWO_PI*val[0];
2792  Double_t a_w0 = a/w0;
2793  Double_t w0t = w0*t;
2794 
2795  Double_t j1, j0;
2796  if (fabs(w0t) < 0.001) { // check zero time limits of the spherical bessel functions j0(x) and j1(x)
2797  j0 = 1.0;
2798  j1 = 0.0;
2799  } else {
2800  j0 = sin(w0t)/w0t;
2801  j1 = (sin(w0t)-w0t*cos(w0t))/(w0t*w0t);
2802  }
2803 
2804  p0exp[i] = 1.0 - a_w0*j1*exp(-at) - a_w0*a_w0*(j0*exp(-at) - 1.0) - GetLFIntegralValue(t);
2805  }
2806  break;
2807  default:
2808  break;
2809  }
2810  p0exp[i] *= TMath::Exp(-val[2]*t);
2811  t += dt;
2812  }
2813 
2814  // solve the volterra equation (trapezoid integration)
2815  fDynLFFuncValue[0]=p0exp[0];
2816 
2817  Double_t sum;
2818  Double_t a;
2819  Double_t preFactor = dt*val[2];
2820  for (UInt_t i=1; i<N; i++) {
2821  sum = p0exp[i];
2822  sum += 0.5*preFactor*p0exp[i]*fDynLFFuncValue[0];
2823  for (UInt_t j=1; j<i; j++) {
2824  sum += preFactor*p0exp[i-j]*fDynLFFuncValue[j];
2825  }
2826  a = 1.0-0.5*preFactor*p0exp[0];
2827 
2828  fDynLFFuncValue[i]=sum/a;
2829  }
2830 
2831  // clean up
2832  p0exp.clear();
2833 }
2834 
2835 //--------------------------------------------------------------------------
2843 Double_t PTheory::GetDynKTLFValue(const Double_t t) const
2844 {
2845  if (t < 0.0)
2846  return 0.0;
2847 
2848  UInt_t idx = static_cast<UInt_t>(t/fDynLFdt);
2849 
2850  if (idx + 2 > fDynLFFuncValue.size())
2851  return fDynLFFuncValue.back();
2852 
2853  // linearly interpolate between the two relvant function bins
2854  Double_t df = (fDynLFFuncValue[idx+1]-fDynLFFuncValue[idx])*(t/fDynLFdt-static_cast<Double_t>(idx));
2855 
2856  return fDynLFFuncValue[idx]+df;
2857 }
2858 
2859 //--------------------------------------------------------------------------
2873 Double_t PTheory::MuMinusExpTF(Double_t t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
2874 {
2875  // expected parameters: N0 tau A lambda phase frequency [tshift]
2876 
2877  Double_t val[7];
2878 
2879  assert(fParamNo.size() <= 7);
2880 
2881  // check if FUNCTIONS are used
2882  for (UInt_t i=0; i<fParamNo.size(); i++) {
2883  if (fParamNo[i] < MSR_PARAM_FUN_OFFSET) { // parameter or resolved map
2884  val[i] = paramValues[fParamNo[i]];
2885  } else { // function
2886  val[i] = funcValues[fParamNo[i]-MSR_PARAM_FUN_OFFSET];
2887  }
2888  }
2889 
2890  Double_t tt;
2891  if (fParamNo.size() == 6) // no tshift
2892  tt = t;
2893  else // tshift present
2894  tt = t-val[6];
2895 
2896  return val[0]*exp(-tt/val[1])*(1.0+val[2]*exp(-val[3]*tt)*cos(TWO_PI*val[5]*tt+DEG_TO_RAD*val[4]));
2897 }
2898 
2899 //--------------------------------------------------------------------------
2900 // END
2901 //--------------------------------------------------------------------------
PTheory::StaticGaussKTLF
virtual Double_t StaticGaussKTLF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1214
PTheory::fUserFcn
PUserFcnBase * fUserFcn
pointer to the user function object
Definition: PTheory.h:300
PTheory::PTheory
PTheory(PMsrHandler *msrInfo, UInt_t runNo, const Bool_t hasParent=false)
Definition: PTheory.cpp:93
PTheory::DynamicNKZF
virtual Double_t DynamicNKZF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2363
THEORY_GENERAL_EXP
#define THEORY_GENERAL_EXP
Definition: PTheory.h:49
THEORY_CONST
#define THEORY_CONST
Definition: PTheory.h:46
PTheory::StaticGaussKT
virtual Double_t StaticGaussKT(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1172
PTheory::StrKT
virtual Double_t StrKT(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1697
THEORY_BESSEL
#define THEORY_BESSEL
Definition: PTheory.h:66
PTheory::fType
UInt_t fType
function tag
Definition: PTheory.h:292
DEG_TO_RAD
#define DEG_TO_RAD
Definition: PTheory.h:115
PTheory::MakeCleanAndTidyTheoryBlock
virtual void MakeCleanAndTidyTheoryBlock(PMsrLines *fullTheoryBlock)
Definition: PTheory.cpp:794
THEORY_INTERNAL_FIELD_KORNILOV
#define THEORY_INTERNAL_FIELD_KORNILOV
Definition: PTheory.h:64
PTheory::Polynom
virtual Double_t Polynom(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2478
PUserFcnBase::GlobalPartIsValid
virtual Bool_t GlobalPartIsValid() const
if a user function is using a global part, this function returns if the global object part is valid (...
Definition: PUserFcnBase.h:50
PTheory::SearchDataBase
virtual Int_t SearchDataBase(TString name)
Definition: PTheory.cpp:739
PTheory::SkewedGauss
virtual Double_t SkewedGauss(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2163
PTheory::MakeCleanAndTidyPolynom
virtual void MakeCleanAndTidyPolynom(UInt_t i, PMsrLines *fullTheoryBlock)
Definition: PTheory.cpp:885
PTheory::Bessel
virtual Double_t Bessel(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2073
PTheory::DynamicLorentzKTLF
virtual Double_t DynamicLorentzKTLF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1541
THEORY_ABRAGAM
#define THEORY_ABRAGAM
Definition: PTheory.h:61
THEORY_TF_COS
#define THEORY_TF_COS
Definition: PTheory.h:62
THEORY_STATIC_GAUSS_KT
#define THEORY_STATIC_GAUSS_KT
Definition: PTheory.h:51
theo_data_base::fComment
TString fComment
comment added in the msr-file theory block to help the used
Definition: PTheory.h:131
THEORY_ASYMMETRY
#define THEORY_ASYMMETRY
Definition: PTheory.h:47
PTheory::fMul
PTheory * fMul
pointers to the add-sub-function or the multiply-sub-function
Definition: PTheory.h:295
PTheory::MuMinusExpTF
virtual Double_t MuMinusExpTF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2873
THEORY_STATIC_LORENTZ_KT
#define THEORY_STATIC_LORENTZ_KT
Definition: PTheory.h:54
PTheory::fUserFcnClassName
TString fUserFcnClassName
name of the user function class for within root
Definition: PTheory.h:298
THEORY_MAX_PARAM
#define THEORY_MAX_PARAM
Definition: PTheory.h:112
THEORY_DYNAMIC_GAUSS_KT_LF
#define THEORY_DYNAMIC_GAUSS_KT_LF
Definition: PTheory.h:53
THEORY_SIMPLE_GAUSS
#define THEORY_SIMPLE_GAUSS
Definition: PTheory.h:50
THEORY_UNDEFINED
#define THEORY_UNDEFINED
Definition: PTheory.h:45
THEORY_MAX
#define THEORY_MAX
Definition: PTheory.h:109
PTheory::fUserParam
PDoubleVector fUserParam
vector holding the resolved user function parameters, i.e. map and function resolved.
Definition: PTheory.h:301
PTheory::InternalField
virtual Double_t InternalField(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1927
THEORY_POLYNOM
#define THEORY_POLYNOM
Definition: PTheory.h:74
theo_data_base::fType
UInt_t fType
function tag
Definition: PTheory.h:126
THEORY_STR_KT
#define THEORY_STR_KT
Definition: PTheory.h:58
PTheory::fNoOfParam
UInt_t fNoOfParam
number of parameters for the given function
Definition: PTheory.h:294
PIntVector
std::vector< Int_t > PIntVector
Definition: PMusr.h:178
SQRT_PI
#define SQRT_PI
Definition: PTheory.cpp:47
PMsrHandler::GetFuncIndex
virtual UInt_t GetFuncIndex(Int_t funNo)
Definition: PMsrHandler.h:95
PTheory::Func
virtual Double_t Func(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:390
PTheory::StaticNKTF
virtual Double_t StaticNKTF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2307
PTheory::UserFcn
virtual Double_t UserFcn(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2515
PDoubleVector
std::vector< Double_t > PDoubleVector
Definition: PMusr.h:196
PTheory::GetUserFcnIdx
virtual Int_t GetUserFcnIdx(UInt_t lineNo) const
Definition: PTheory.cpp:765
THEORY_STATIC_TF_NK
#define THEORY_STATIC_TF_NK
Definition: PTheory.h:70
PTheory::TFCos
virtual Double_t TFCos(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1886
PMsrHandler::GetMsrRunList
virtual PMsrRunList * GetMsrRunList()
Definition: PMsrHandler.h:63
MSR_PARAM_FUN_OFFSET
#define MSR_PARAM_FUN_OFFSET
Definition: PMusr.h:127
theo_data_base::fName
TString fName
name of the function as written into the msr-file
Definition: PTheory.h:129
PTheory::CalculateGaussLFIntegral
virtual void CalculateGaussLFIntegral(const Double_t *val) const
Definition: PTheory.cpp:2542
PTheory::CalculateLorentzLFIntegral
virtual void CalculateLorentzLFIntegral(const Double_t *val) const
Definition: PTheory.cpp:2612
PTheory::fUserFcnIdx
Int_t fUserFcnIdx
index of the user function within the theory tree
Definition: PTheory.h:297
PTheory::fAdd
PTheory * fAdd
Definition: PTheory.h:295
PUserFcnBase::NeedGlobalPart
virtual Bool_t NeedGlobalPart() const
if a user function needs a global part this function should return true, otherwise false (default: fa...
Definition: PUserFcnBase.h:48
THEORY_STATIC_ZF_NK
#define THEORY_STATIC_ZF_NK
Definition: PTheory.h:69
PTheory::fUserFcnSharedLibName
TString fUserFcnSharedLibName
name of the shared lib to which the user function belongs
Definition: PTheory.h:299
PTheory::GetDynKTLFValue
virtual Double_t GetDynKTLFValue(const Double_t t) const
Definition: PTheory.cpp:2843
THEORY_DYNAMIC_ZF_NK
#define THEORY_DYNAMIC_ZF_NK
Definition: PTheory.h:71
PTheory::SimpleGauss
virtual Double_t SimpleGauss(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1132
PTheory::fDynLFFuncValue
PDoubleVector fDynLFFuncValue
needed for LF. Keeps the dynamic LF KT function values
Definition: PTheory.h:309
PTheory::DynamicGaussKTLF
virtual Double_t DynamicGaussKTLF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1303
THEORY_MU_MINUS_EXP
#define THEORY_MU_MINUS_EXP
Definition: PTheory.h:73
PMsrLineStructure::fLineNo
Int_t fLineNo
original line number of the msr-file
Definition: PMusr.h:531
THEORY_RANDOM_ANISOTROPIC_HYPERFINE
#define THEORY_RANDOM_ANISOTROPIC_HYPERFINE
Definition: PTheory.h:60
PTheory::StaticLorentzKTLF
virtual Double_t StaticLorentzKTLF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1444
THEORY_USER_FCN
#define THEORY_USER_FCN
Definition: PTheory.h:75
THEORY_COMBI_LGKT
#define THEORY_COMBI_LGKT
Definition: PTheory.h:57
PMsrLines
std::vector< PMsrLineStructure > PMsrLines
Definition: PMusr.h:539
PTheory::GetLFIntegralValue
virtual Double_t GetLFIntegralValue(const Double_t t) const
Definition: PTheory.cpp:2682
PMsrLineStructure::fLine
TString fLine
msr-file line
Definition: PMusr.h:532
PTheory::StaticNKZF
virtual Double_t StaticNKZF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2252
gGlobalUserFcn
std::vector< void * > gGlobalUserFcn
Definition: PUserFcnBase.cpp:80
PTheory::MakeCleanAndTidyUserFcn
virtual void MakeCleanAndTidyUserFcn(UInt_t i, PMsrLines *fullTheoryBlock)
Definition: PTheory.cpp:942
PTheory::fMsrInfo
PMsrHandler * fMsrInfo
pointer to the msr-file handler
Definition: PTheory.h:303
PTheory::InternalFieldLL
virtual Double_t InternalFieldLL(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2021
PUserFcnBase::SetGlobalPart
virtual void SetGlobalPart(std::vector< void *> &globalPart, UInt_t idx)
if a user function is using a global part, this function is used to invoke and retrieve the proper gl...
Definition: PUserFcnBase.h:49
THEORY_STATIC_LORENTZ_KT_LF
#define THEORY_STATIC_LORENTZ_KT_LF
Definition: PTheory.h:55
PTheory::~PTheory
virtual ~PTheory()
Definition: PTheory.cpp:332
PTheory::IsValid
virtual Bool_t IsValid()
Definition: PTheory.cpp:361
PTheory::GeneralExp
virtual Double_t GeneralExp(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1084
PTheory::SimpleExp
virtual Double_t SimpleExp(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1044
PTheory::CalculateDynKTLF
virtual void CalculateDynKTLF(const Double_t *val, Int_t tag) const
Definition: PTheory.cpp:2709
PTheory::fSamplingTime
Double_t fSamplingTime
needed for LF. Keeps the sampling time of the non-analytic integral
Definition: PTheory.h:305
PTheory::DynamicNKTF
virtual Double_t DynamicNKTF(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2424
PTheory::Asymmetry
virtual Double_t Asymmetry(const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1014
PMsrHandler::GetMsrTheory
virtual PMsrLines * GetMsrTheory()
Definition: PMsrHandler.h:60
SQRT_TWO
#define SQRT_TWO
Definition: PTheory.cpp:46
THEORY_DYNAMIC_LORENTZ_KT_LF
#define THEORY_DYNAMIC_LORENTZ_KT_LF
Definition: PTheory.h:56
PTheory::Abragam
virtual Double_t Abragam(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1842
PTheory::InternalBessel
virtual Double_t InternalBessel(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:2114
PTheory::fValid
Bool_t fValid
flag to tell if the theory is valid
Definition: PTheory.h:291
PTheory::RandomAnisotropicHyperfine
virtual Double_t RandomAnisotropicHyperfine(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1797
PTheory::CombiLGKT
virtual Double_t CombiLGKT(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1651
PTheory::InternalFieldGK
virtual Double_t InternalFieldGK(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1969
fgTheoDataBase
static PTheoDataBase fgTheoDataBase[THEORY_MAX]
Definition: PTheory.h:139
THEORY_STATIC_GAUSS_KT_LF
#define THEORY_STATIC_GAUSS_KT_LF
Definition: PTheory.h:52
THEORY_INTERNAL_FIELD
#define THEORY_INTERNAL_FIELD
Definition: PTheory.h:63
theo_data_base::fNoOfParam
UInt_t fNoOfParam
number of parameters for this function
Definition: PTheory.h:127
PTheory::fLFIntegral
PDoubleVector fLFIntegral
needed for LF. Keeps the non-analytic integral values
Definition: PTheory.h:307
THEORY_SPIN_GLASS
#define THEORY_SPIN_GLASS
Definition: PTheory.h:59
THEORY_SKEWED_GAUSS
#define THEORY_SKEWED_GAUSS
Definition: PTheory.h:68
status
return status
Definition: PStartupHandler.cpp:85
THEORY_INTERNAL_FIELD_LARKIN
#define THEORY_INTERNAL_FIELD_LARKIN
Definition: PTheory.h:65
TWO_PI
#define TWO_PI
Definition: PTheory.h:117
THEORY_INTERNAL_BESSEL
#define THEORY_INTERNAL_BESSEL
Definition: PTheory.h:67
PTheory::fParamNo
std::vector< UInt_t > fParamNo
holds the parameter numbers for the theory (including maps and functions, see constructor desciption)...
Definition: PTheory.h:293
PTheory::fPrevParam
Double_t fPrevParam[THEORY_MAX_PARAM]
needed for LF. Keeps the previous fitting parameters
Definition: PTheory.h:306
PTheory::SpinGlass
virtual Double_t SpinGlass(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1746
PTheory::Constant
virtual Double_t Constant(const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:986
THEORY_DYNAMIC_TF_NK
#define THEORY_DYNAMIC_TF_NK
Definition: PTheory.h:72
THEORY_SIMPLE_EXP
#define THEORY_SIMPLE_EXP
Definition: PTheory.h:48
PTheory::CleanUp
virtual void CleanUp(PTheory *theo)
Definition: PTheory.cpp:716
PTheory::fDynLFdt
Double_t fDynLFdt
needed for LF. Keeps the time step for the dynamic LF calculation, used in the integral equation appr...
Definition: PTheory.h:308
PTheory::StaticLorentzKT
virtual Double_t StaticLorentzKT(Double_t t, const PDoubleVector &paramValues, const PDoubleVector &funcValues) const
Definition: PTheory.cpp:1401
theo_data_base::fCommentTimeShift
TString fCommentTimeShift
comment added in the msr-file theory block if there is a time shift
Definition: PTheory.h:132