Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
TFormula_v5.cxx
Go to the documentation of this file.
1 // @(#)root/hist:$Id$
2 // Author: Nicolas Brun 19/08/95
3 
4 /*************************************************************************
5  * Copyright (C) 1995-2000, Rene Brun and Fons Rademakers. *
6  * All rights reserved. *
7  * *
8  * For the licensing terms see $ROOTSYS/LICENSE. *
9  * For the list of contributors see $ROOTSYS/README/CREDITS. *
10  *************************************************************************/
11 
12 #include <cmath>
13 
14 #include "Riostream.h"
15 #include "TROOT.h"
16 #include "TClass.h"
17 #include "v5/TFormula.h"
18 #include "TMath.h"
19 #include "TRandom.h"
20 #include "TFunction.h"
21 #include "TMethodCall.h"
22 #include "TObjString.h"
23 #include "TError.h"
24 #include "v5/TFormulaPrimitive.h"
25 #include "TInterpreter.h"
26 #include "TVirtualMutex.h"
27 
28 #ifdef WIN32
29 #pragma optimize("",off)
30 #endif
31 
32 static Int_t gMAXOP=1000,gMAXPAR=1000,gMAXCONST=1000;
33 const Int_t gMAXSTRINGFOUND = 10;
34 const UInt_t kOptimizationError = BIT(19);
35 
36 ClassImp(ROOT::v5::TFormula);
37 
38 namespace ROOT {
39 
40  namespace v5 {
41 
42 /** \class TFormula TFormula.h "inc/v5/TFormula.h"
43  \ingroup Hist
44 The FORMULA class (ROOT version 5)
45 
46  Example of valid expressions:
47 
48  - `sin(x)/x`
49  - `[0]*sin(x) + [1]*exp(-[2]*x)`
50  - `x + y**2`
51  - `x^2 + y^2`
52  - `[0]*pow([1],4)`
53  - `2*pi*sqrt(x/y)`
54  - `gaus(0)*expo(3) + ypol3(5)*x`
55  - `gausn(0)*expo(3) + ypol3(5)*x`
56 
57  In the last example above:
58 
59  gaus(0) is a substitute for `[0]*exp(-0.5*((x-[1])/[2])**2)`
60  and (0) means start numbering parameters at 0
61 
62  gausn(0) is a substitute for `[0]*exp(-0.5*((x-[1])/[2])**2)/(sqrt(2*pi)*[2]))`
63  and (0) means start numbering parameters at 0
64 
65  expo(3) is a substitute for `exp([3]+[4]*x)`
66 
67  pol3(5) is a substitute for `par[5]+par[6]*x+par[7]*x**2+par[8]*x**3`
68  (here Pol3 stands for Polynomial of degree 3)
69 
70  TMath functions can be part of the expression, eg:
71 
72  - `TMath::Landau(x)*sin(x)`
73  - `TMath::Erf(x)`
74 
75  Comparisons operators are also supported (&&, ||, ==, <=, >=, !)
76  Examples:
77 
78  sin(x*(x<0.5 || x>1))
79 
80  If the result of a comparison is TRUE, the result is 1, otherwise 0.
81 
82  Already predefined names can be given. For example, if the formula
83 
84  TFormula old(sin(x*(x<0.5 || x>1))) one can assign a name to the formula. By default
85  the name of the object = title = formula itself.
86  old.SetName("old").
87  then, old can be reused in a new expression.
88  TFormula new("x*old") is equivalent to:
89  TFormula new("x*sin(x*(x<0.5 || x>1))")
90 
91  Up to 4 dimensions are supported (indicated by x, y, z, t)
92  An expression may have 0 parameters or a list of parameters
93  indicated by the sequence [par_number]
94 
95  A graph showing the logic to compile and analyze a formula
96  is shown in TFormula::Compile and TFormula::Analyze.
97  Once a formula has been compiled, it can be evaluated for a given
98  set of parameters. see graph in TFormula::EvalPar.
99 
100  This class is the base class for the function classes TF1,TF2 and TF3.
101  It is also used by the ntuple selection mechanism TNtupleFormula.
102 
103  In version 7 of TFormula, the usage of fOper has been changed
104  to improve the performance of TFormula::EvalPar.
105  Conceptually, fOper was changed from a simple array of Int_t
106  to an array of composite values.
107  For example a 'ylandau(5)' operation used to be encoded as 4105;
108  it is now encoded as (klandau >> kTFOperShift) + 5
109  Any class inheriting from TFormula and using directly fOper (which
110  is now a private data member), needs to be updated to take this
111  in consideration. The member functions recommended to set and
112  access fOper are: SetAction, GetAction, GetActionParam
113  For more performant access to the information, see the implementation
114  TFormula::EvalPar
115 
116  ### CHANGING DEFAULT SETTINGS
117 
118  When creating complex formula , it may be necessary to increase
119  some default parameters. see static function TFormula::SetMaxima
120 
121  ### WHY TFormula CANNOT ACCEPT A CLASS MEMBER FUNCTION ?
122 
123  This is a frequently asked question.
124  C++ is a strongly typed language. There is no way for TFormula (without
125  recompiling this class) to know about all possible user defined data types.
126  This also apply to the case of a static class function.
127  Because TMath is a special and frequent case, TFormula is aware
128  of all TMath functions.
129 */
130 
131 ////////////////////////////////////////////////////////////////////////////////
132 /// Formula default constructor.
133 
134 TFormula::TFormula(): TNamed()
135 {
136  fNdim = 0;
137  fNpar = 0;
138  fNoper = 0;
139  fNconst = 0;
140  fNumber = 0;
141  fExpr = 0;
142  fOper = 0;
143  fConst = 0;
144  fParams = 0;
145  fNstring= 0;
146  fNames = 0;
147  fNval = 0;
148  //
149  //MI change
150  fNOperOptimized = 0;
151  fExprOptimized = 0;
152  fOperOptimized = 0;
153  fOperOffset = 0;
154  fPredefined = 0;
155  fOptimal = (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalParOld;
156 }
157 
158 ////////////////////////////////////////////////////////////////////////////////
159 /// Normal Formula constructor.
160 
161 TFormula::TFormula(const char *name,const char *expression) :
162  TNamed(name,expression)
163 {
164  fNdim = 0;
165  fNpar = 0;
166  fNoper = 0;
167  fNconst = 0;
168  fNumber = 0;
169  fExpr = 0;
170  fOper = 0;
171  fConst = 0;
172  fParams = 0;
173  fNstring= 0;
174  fNames = 0;
175  fNval = 0;
176  //
177  //MI change
178  fNOperOptimized = 0;
179  fExprOptimized = 0;
180  fOperOptimized = 0;
181  fOperOffset = 0;
182  fPredefined = 0;
183  fOptimal = (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalParOld;
184 
185  if (!expression || !*expression) {
186  Error("TFormula", "expression may not be 0 or have 0 length");
187  return;
188  }
189 
190  //eliminate blanks in expression
191  Int_t i,j,nch;
192  nch = strlen(expression);
193  char *expr = new char[nch+1];
194  j = 0;
195  for (i=0;i<nch;i++) {
196  if (expression[i] == ' ') continue;
197  if (i > 0 && (expression[i] == '*') && (expression[i-1] == '*')) {
198  expr[j-1] = '^';
199  continue;
200  }
201  expr[j] = expression[i]; j++;
202  }
203  expr[j] = 0;
204  Bool_t gausNorm = kFALSE;
205  Bool_t landauNorm = kFALSE;
206  Bool_t linear = kFALSE;
207 
208  if (j) {
209  TString chaine = expr;
210  //special case for functions for linear fitting
211  if (chaine.Contains("++"))
212  linear = kTRUE;
213  // special case for normalized gaus
214  if (chaine.Contains("gausn")) {
215  gausNorm = kTRUE;
216  TString tmp = chaine;
217  tmp.ReplaceAll("gausn","");
218  tmp.ReplaceAll("landaun","");
219  if ( tmp.Contains("gaus") )
220  Warning("TFormula","Cannot use both gaus and gausn - gaus will be treated as gausn");
221  if ( tmp.Contains("landau") )
222  Warning("TFormula","Cannot use both gausn and landau - landau will be treated as landaun");
223  }
224  // special case for normalized landau
225  if (chaine.Contains("landaun")) {
226  landauNorm = kTRUE;
227  TString tmp = chaine;
228  tmp.ReplaceAll("landaun","");
229  tmp.ReplaceAll("gausn","");
230  if ( tmp.Contains("gaus") ) {
231  Warning("TFormula","Cannot use both gaus and landaun - gaus will be treated as gausn");
232  }
233  if ( tmp.Contains("landau") )
234  Warning("TFormula","Cannot use both landau and landaun - landau will be treated as landaun");
235  }
236  // need to the replacement here for the error message before
237  if (gausNorm)
238  chaine.ReplaceAll("gausn","gaus");
239  if (landauNorm)
240  chaine.ReplaceAll("landaun","landau");
241 
242  SetTitle(chaine.Data());
243  }
244  delete [] expr;
245 
246  if (linear) SetBit(kLinear);
247 
248  if (Compile()) return;
249 
250  if (gausNorm) SetBit(kNormalized);
251  if (landauNorm) SetBit(kNormalized);
252 
253  // Store formula in linked list of formula in ROOT
254 
255 
256  if (strcmp(name,"x")==0 || strcmp(name,"y")==0 ||
257  strcmp(name,"z")==0 || strcmp(name,"t")==0 )
258  {
259  Error("TFormula","The name \'%s\' is reserved as a TFormula variable name.\n"
260  "\tThis function will not be registered in the list of functions",name);
261  } else {
262  R__LOCKGUARD(gROOTMutex);
263  TFormula *old = (TFormula*)gROOT->GetListOfFunctions()->FindObject(name);
264  if (old) {
265  gROOT->GetListOfFunctions()->Remove(old);
266  }
267  gROOT->GetListOfFunctions()->Add(this);
268  }
269 }
270 
271 ////////////////////////////////////////////////////////////////////////////////
272 /// Default constructor.
273 
274 TFormula::TFormula(const TFormula &formula) : TNamed()
275 {
276  fNdim = 0;
277  fNpar = 0;
278  fNoper = 0;
279  fNconst = 0;
280  fNumber = 0;
281  fExpr = 0;
282  fOper = 0;
283  fConst = 0;
284  fParams = 0;
285  fNstring= 0;
286  fNames = 0;
287  fNval = 0;
288  fNOperOptimized = 0;
289  fPredefined = 0;
290  fOperOffset = 0;
291  fExprOptimized = 0;
292  fOperOptimized = 0;
293  fOptimal = (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalParOld;
294 
295  ((TFormula&)formula).TFormula::Copy(*this);
296 }
297 
298 ////////////////////////////////////////////////////////////////////////////////
299 /// Operator =
300 
301 TFormula& TFormula::operator=(const TFormula &rhs)
302 {
303  if (this != &rhs) {
304  rhs.Copy(*this);
305  }
306  return *this;
307 }
308 
309 ////////////////////////////////////////////////////////////////////////////////
310 /// Formula default destructor.
311 
312 TFormula::~TFormula()
313 {
314  if (gROOT) {
315  R__LOCKGUARD(gROOTMutex);
316  gROOT->GetListOfFunctions()->Remove(this);
317  }
318 
319  ClearFormula();
320 }
321 
322 ////////////////////////////////////////////////////////////////////////////////
323 /// Check if the chain as function call.
324 ///
325 /// If you overload this member function, you also HAVE TO
326 /// never call the constructor:
327 ///
328 /// ~~~ {.cpp}
329 /// TFormula::TFormula(const char *name,const char *expression)
330 /// ~~~
331 ///
332 /// and write your own constructor
333 ///
334 /// ~~~ {.cpp}
335 /// MyClass::MyClass(const char *name,const char *expression) : TFormula()
336 /// ~~~
337 ///
338 /// which has to call the TFormula default constructor and whose implementation
339 /// should be similar to the implementation of the normal TFormula constructor
340 ///
341 /// This is necessary because the normal TFormula constructor call indirectly
342 /// the virtual member functions Analyze, DefaultString, DefaultValue
343 /// and DefaultVariable.
344 
345 Bool_t TFormula::AnalyzeFunction(TString &chaine, Int_t &err, Int_t offset)
346 {
347  int i;
348 
349  // We have to decompose the chain is 3 potential components:
350  // namespace::functionName( args )
351 
352  Ssiz_t argStart = chaine.First('(');
353  if (argStart<0) return false;
354 
355  TString functionName = chaine(0,argStart);
356 
357  // This does not support template yet (where the scope operator might be in
358  // one of the template arguments
359  Ssiz_t scopeEnd = functionName.Last(':');
360  TString spaceName;
361  if (scopeEnd>0 && functionName[scopeEnd-1]==':') {
362  spaceName = functionName(0,scopeEnd-1);
363  functionName.Remove(0,scopeEnd+1);
364  }
365 
366  // Now we need to count and decompose the actual arguments, we could also check the type
367  // of the arguments
368  if (chaine[chaine.Length()-1] != ')') {
369  Error("AnalyzeFunction","We thought we had a function but we dont (in %s)\n",chaine.Data());
370  }
371 
372  TString args = chaine(argStart+1,chaine.Length()-2-argStart);
373  TObjArray argArr;
374  argArr.SetOwner(kTRUE);
375  //fprintf(stderr,"args are '%s'\n",args.Data());
376 
377  Bool_t inString = false;
378  int paran = 0;
379  int brack = 0;
380  int prevComma = 0;
381  int nargs = 0;
382  for(i=0; i<args.Length(); i++) {
383  if (args[i]=='"') inString = !inString;
384  if (inString) continue;
385 
386  Bool_t foundArg = false;
387  switch(args[i]) {
388 
389  case '(': paran++; break;
390  case ')': paran--; break;
391  case '[': brack++; break;
392  case ']': brack--; break;
393 
394  case ',': if (paran==0 && brack==0) { foundArg = true; } break;
395  }
396  if ((i+1)==args.Length()) {
397  foundArg = true; i++;
398  }
399  if (foundArg) {
400  TString arg = args(prevComma,i-prevComma);
401 
402  // Here we could
403  // a) check the type
404  //fprintf(stderr,"found #%d arg %s\n",nargs,arg.Data());
405 
406  // We register the arg for later usage
407  argArr.Add(new TObjString(arg));
408  nargs++;
409 
410  prevComma = i+1;
411  };
412  }
413 
414  if (nargs>999) {
415  err = 7;
416  return false;
417  }
418 
419  // Now we need to lookup the function and check its arguments.
420  TClass *ns = (spaceName.Length()) ? TClass::GetClass(spaceName) : 0;
421  ClassInfo_t *cinfo = 0;
422  if (ns) {
423  cinfo = ns->GetClassInfo();
424  } else {
425  cinfo = gInterpreter->ClassInfo_Factory();
426  }
427 
428  // ROOT does yet have a complete TType class, but TCling does,
429  // so let's use that for now.
430  static TypeInfo_t *const doubletype { gInterpreter->TypeInfo_Factory("double") };
431 
432  std::vector<TypeInfo_t*> proto(nargs,doubletype);
433 
434  CallFunc_t *callfunc = gInterpreter->CallFunc_Factory();
435  Long_t func_offset;
436  gInterpreter->CallFunc_SetFuncProto(callfunc,cinfo,functionName,proto,false,&func_offset,ROOT::kConversionMatch);
437 
438  TMethodCall *method = new TMethodCall(ns,callfunc,func_offset);
439 
440  if (!ns) gInterpreter->ClassInfo_Delete(cinfo);
441  gInterpreter->CallFunc_Delete(callfunc);
442 
443  if (method->IsValid()) {
444  if (method->ReturnType() == TMethodCall::kOther) {
445  /*
446  Error("Compile",
447  "TFormula can only call interpreted and compiled function that returns a numerical type %s returns a %s\n",
448  method->GetMethodName(), method->GetMethod()->GetReturnTypeName());
449  */
450  err=29;
451 
452  } else {
453 
454  // Analyze the arguments
455  TIter next(&argArr);
456  TObjString *objstr;
457  while ( (objstr=(TObjString*)next()) ) {
458  Analyze(objstr->String(),err,offset);
459  }
460 
461  fFunctions.Add(method);
462  fExpr[fNoper] = method->GetMethod()->GetPrototype();
463  SetAction(fNoper, kFunctionCall, fFunctions.GetLast()*1000 + nargs);
464  fNoper++;
465  return true;
466  }
467  }
468 
469  delete method;
470  //
471  // MI change - extended space of functions
472  // not forward compatible change
473  //
474  TString cbase(chaine);
475  Int_t args_paran = cbase.First("(");
476  if (args_paran>0){
477  cbase[args_paran]=0;
478  }
479 
480  ROOT::v5::TFormulaPrimitive *prim = ROOT::v5::TFormulaPrimitive::FindFormula(cbase, args_paran>0 ? cbase.Data() + args_paran + 1 : (const char*)0);
481  if (prim && (!IsA()->GetBaseClass("TTreeFormula"))) {
482  // TO BE DONE ALSO IN TTREFORMULA - temporary fix MI
483  // Analyze the arguments
484  TIter next(&argArr);
485  TObjString *objstr;
486  while ( (objstr=(TObjString*)next()) ) {
487  Analyze(objstr->String(),err,offset); if (err) return kFALSE;
488  }
489  if (nargs!=prim->fNArguments) {
490  Error("Compile", "%s requires %d arguments",
491  prim->GetName(), prim->fNArguments);
492  return kFALSE;
493  }
494  fExpr[fNoper] = prim->GetName();
495  if (prim->fType==10){
496  SetAction(fNoper, kFD1);
497  }
498  if (prim->fType==110){
499  SetAction(fNoper, kFD2);
500  }
501  if (prim->fType==1110){
502  SetAction(fNoper, kFD3);
503  }
504  if (prim->fType==-1){
505  SetAction(fNoper, kFDM);
506  if (fNpar<prim->fNParameters) fNpar+=prim->fNParameters;
507  }
508 
509  fNoper++;
510  return kTRUE;
511  }
512 
513  return kFALSE;
514 }
515 
516 ////////////////////////////////////////////////////////////////////////////////
517 /// Analyze a sub-expression in one formula.
518 ///
519 /// Expressions in one formula are recursively analyzed.
520 /// Result of analysis is stored in the object tables.
521 ///
522 /// ### Table of function codes and errors
523 ///
524 /// ~~~ {.cpp}
525 /// * functions :
526 ///
527 /// + 1 pow 20
528 /// - 2 sq 21
529 /// * 3 sqrt 22
530 /// / 4 strstr 23
531 /// % 5 min 24
532 /// max 25
533 /// log 30
534 /// cos 10 exp 31
535 /// sin 11 log10 32
536 /// tan 12
537 /// acos 13 abs 41
538 /// asin 14 sign 42
539 /// atan 15 int 43
540 /// atan2 16
541 /// fmod 17 rndm 50
542 ///
543 /// cosh 70 acosh 73
544 /// sinh 71 asinh 74
545 /// tanh 72 atanh 75
546 ///
547 /// expo 100 gaus 110 gausn (see note below)
548 /// expo(0) 100 0 gaus(0) 110 0 gausn(0)
549 /// expo(1) 100 1 gaus(1) 110 1 gausn(1)
550 /// xexpo 100 x xgaus 110 x xgausn
551 /// yexpo 101 x ygaus 111 x ygausn
552 /// zexpo 102 x zgaus 112 x zgausn
553 /// xyexpo 105 x xygaus 115 x xygausn
554 /// yexpo(5) 102 5 ygaus(5) 111 5 ygausn(5)
555 /// xyexpo(2) 105 2 xygaus(2) 115 2 xygausn(2)
556 ///
557 /// landau 120 x landaun (see note below)
558 /// landau(0) 120 0 landaun(0)
559 /// landau(1) 120 1 landaun(1)
560 /// xlandau 120 x xlandaun
561 /// ylandau 121 x ylandaun
562 /// zlandau 122 x zlandaun
563 /// xylandau 125 x xylandaun
564 /// ylandau(5) 121 5 ylandaun(5)
565 /// xylandau(2) 125 2 xylandaun(2)
566 ///
567 /// pol0 130 x pol1 130 1xx
568 /// pol0(0) 130 0 pol1(0) 130 100
569 /// pol0(1) 130 1 pol1(1) 130 101
570 /// xpol0 130 x xpol1 130 101
571 /// ypol0 131 x ypol1 131 101
572 /// zpol0 132 x zpol1 132 1xx
573 /// ypol0(5) 131 5 ypol1(5) 131 105
574 ///
575 /// pi 40
576 ///
577 /// && 60 < 64
578 /// || 61 > 65
579 /// == 62 <= 66
580 /// != 63 => 67
581 /// ! 68
582 /// ==(string) 76 & 78
583 /// !=(string) 77 | 79
584 /// <<(shift) 80 >>(shift) 81
585 /// ? : 82
586 ///
587 /// * constants (kConstants) :
588 ///
589 /// c0 141 1 c1 141 2 etc..
590 ///
591 /// * strings (kStringConst):
592 ///
593 /// sX 143 x
594 ///
595 /// * variables (kFormulaVar) :
596 ///
597 /// x 144 0 y 144 1 z 144 2 t 144 3
598 ///
599 /// * parameters :
600 ///
601 /// [1] 140 1
602 /// [2] 140 2
603 /// etc.
604 /// ~~~
605 ///
606 /// ### Special cases for normalized gaussian or landau distributions
607 ///
608 /// the expression "gaus" is a substitute for
609 ///
610 /// [0]*exp(-0.5*((x-[1])/[2])**2)
611 ///
612 /// to obtain a standard normalized gaussian, use "gausn" instead of "gaus"
613 /// the expression "gausn" is a substitute for
614 ///
615 /// [0]*exp(-0.5*((x-[1])/[2])**2)/(sqrt(2*pi)*[2]))
616 ///
617 /// WARNING: gaus and gausn are mutually exclusive in the same expression.
618 ///
619 /// In the same way the expression "landau" is a substitute for
620 ///
621 /// [0]*TMath::Landau(x,[1],[2],kFALSE)
622 ///
623 /// to obtain a standard normalized landau, use "landaun" instead of "landau"
624 /// the expression "landaun" is a substitute for
625 ///
626 /// [0]*TMath::Landau(x,[1],[2],kTRUE)
627 ///
628 /// WARNING: landau and landaun are mutually exclusive in the same expression.
629 ///
630 /// ### Boolean optimization (kBoolOptmize) :
631 ///
632 /// Those pseudo operation are used to implement lazy evaluation of
633 /// && and ||. When the left hand of the expression if false
634 /// (respectively true), the evaluation of the right is entirely skipped
635 /// (since it would not change the value of the expression).
636 ///
637 /// && 142 11 (one operation on right) 142 21 (2 operations on right)
638 /// || 142 12 (one operation on right) 142 22 (2 operations on right)
639 ///
640 /// * functions calls (kFunctionCall) :
641 ///
642 /// f0 145 0 f1 145 1 etc..
643 ///
644 /// ### Errors :
645 ///
646 /// 1 : Division By Zero
647 /// 2 : Invalid Floating Point Operation
648 /// 4 : Empty String
649 /// 5 : invalid syntax
650 /// 6 : Too many operators
651 /// 7 : Too many parameters
652 /// 10 : z specified but not x and y
653 /// 11 : z and y specified but not x
654 /// 12 : y specified but not x
655 /// 13 : z and x specified but not y
656 /// 20 : non integer value for parameter number
657 /// 21 : atan2 requires two arguments
658 /// 22 : pow requires two arguments
659 /// 23 : degree of polynomial not specified
660 /// 24 : Degree of polynomial must be positive
661 /// 25 : Degree of polynomial must be less than 20
662 /// 26 : Unknown name
663 /// 27 : Too many constants in expression
664 /// 28 : strstr requires two arguments
665 /// 29 : interpreted or compiled function have to return a numerical type
666 /// 30 : Bad numerical expression
667 /// 31 : Part of the variable exist but some of it is not accessible or useable
668 /// 40 : '(' is expected
669 /// 41 : ')' is expected
670 /// 42 : '[' is expected
671 /// 43 : ']' is expected
672 ///
673 /// \image html TFormula_analyze.png
674 ///
675 /// ### Special functions
676 ///
677 /// By default, the formula is assigned fNumber=0. However, the following
678 /// formula built with simple functions are assigned fNumber:
679 ///
680 /// "gaus" 100 (or gausn)
681 /// "xygaus" 110
682 /// "expo" 200
683 /// "polN" 300+N
684 /// "landau" 400
685 /// "xylandau" 410
686 ///
687 /// Note that expressions like gaus(0), expo(1) will force fNumber=0
688 ///
689 /// ### Warning when deriving a class from TFormula
690 ///
691 /// If you overload this member function, you also HAVE TO
692 /// never call the constructor:
693 ///
694 /// ~~~ {.cpp}
695 /// TFormula::TFormula(const char *name,const char *expression)
696 /// ~~~
697 ///
698 /// and write your own constructor
699 ///
700 /// ~~~ {.cpp}
701 /// MyClass::MyClass(const char *name,const char *expression) : TFormula()
702 /// ~~~
703 ///
704 /// which has to call the TFormula default constructor and whose implementation
705 /// should be similar to the implementation of the normal TFormula constructor
706 ///
707 /// This is necessary because the normal TFormula constructor call indirectly
708 /// the virtual member functions Analyze, DefaultString, DefaultValue
709 /// and DefaultVariable.
710 
711 void TFormula::Analyze(const char *schain, Int_t &err, Int_t offset)
712 {
713 
714 
715  Int_t valeur,find,n,i,j,k,lchain,nomb,virgule,inter,nest;
716  valeur=find=n=i=j=k=lchain=nomb=virgule=inter=nest = 0;
717  Int_t compt,compt2,compt3,compt4;
718  Bool_t inString;
719  Double_t vafConst;
720  ULong_t vafConst2;
721  Bool_t parenthese;
722  TString s,chaine_error,chaine1ST;
723  TString s1,s2,s3,ctemp;
724 
725  TString chaine = schain;
726  const TFormula *oldformula;
727  Int_t modulo,plus,puiss10,puiss10bis,moins,multi,divi,puiss,et,ou,petit,grand,egal,diff,peteg,grdeg,etx,oux,rshift,lshift,tercond,terelse;
728  char t;
729  TString slash("/"), escapedSlash("\\/");
730  Int_t inter2 = 0;
731  SetNumber(0);
732  Int_t actionCode,actionParam;
733  Int_t err_hint = 0;
734 
735  // Verify correct matching of parenthesis and remove unnecessary parenthesis.
736  lchain = chaine.Length();
737  //if (chaine(lchain-2,2) == "^2") chaine = "sq(" + chaine(0,lchain-2) + ")";
738  parenthese = kTRUE;
739  lchain = chaine.Length();
740  while (parenthese && lchain>0 && err==0){
741  compt = 0;
742  compt2 = 0;
743  inString = false;
744  lchain = chaine.Length();
745  if (lchain==0) err=4;
746  else {
747  for (i=1; i<=lchain; ++i) {
748  if (chaine(i-1,1) == "\"") inString = !inString;
749  if (!inString) {
750  if (chaine(i-1,1) == "[") compt2++;
751  if (chaine(i-1,1) == "]") compt2--;
752  if (chaine(i-1,1) == "(") compt++;
753  if (chaine(i-1,1) == ")") compt--;
754  }
755  if (compt < 0) err = 40; // more open parentheses than close parentheses
756  if (compt2< 0) err = 42; // more ] than [
757  if (compt==0 && (i!=lchain || lchain==1)) parenthese = kFALSE;
758  // if (lchain<3 && chaine(0,1)!="(" && chaine(lchain-1,1)!=")") parenthese = kFALSE;
759  }
760  if (compt > 0) err = 41; // more ( than )
761  if (compt2> 0) err = 43; // more [ than ]
762  if (parenthese) chaine = chaine(1,lchain-2);
763  }
764  } // while parenthesis
765 
766  if (lchain==0) err=4; // empty string
767  modulo=plus=moins=multi=divi=puiss=et=ou=petit=grand=egal=diff=peteg=grdeg=etx=oux=rshift=lshift=tercond=terelse=0;
768 
769  // Look for simple operators
770 
771  if (err==0) {
772  compt = compt2 = compt3 = compt4 = 0;puiss10=0;puiss10bis = 0;
773  inString = false;
774  j = lchain;
775  Bool_t isdecimal = 1; // indicates whether the left part is decimal.
776 
777  for (i=1;i<=lchain; i++) {
778 
779  puiss10=puiss10bis=0;
780  if (i>2) {
781  t = chaine[i-3];
782  isdecimal = isdecimal && (strchr("0123456789.",t)!=0);
783  if (isdecimal) {
784  if ( chaine[i-2] == 'e' || chaine[i-2] == 'E' ) puiss10 = 1;
785  } else if ( strchr("+-/[]()&|><=!*/%^\\",t) ) {
786  isdecimal = 1; // reset after delimiter
787  }
788  }
789  if (j>2) {
790  if (chaine[j-2] == 'e' || chaine[j-2] == 'E') {
791  Bool_t isrightdecimal = 1;
792 
793  for(k=j-3; k>=0 && isrightdecimal; --k) {
794  t = chaine[k];
795  isrightdecimal = isrightdecimal && (strchr("0123456789.",t)!=0);
796  if (!isrightdecimal) {
797  if (strchr("+-/[]()&|><=!*/%^\\",t)!=0) {
798  puiss10bis = 1;
799  }
800  }
801  }
802  if (k<0 && isrightdecimal) puiss10bis = 1;
803  }
804  }
805  if (puiss10 && (i<=lchain)) {
806  t = chaine[i];
807  puiss10 = (strchr("0123456789.",t)!=0);
808  }
809  if (puiss10bis && (j<=lchain)) {
810  t = chaine[j];
811  puiss10bis = (strchr("0123456789.",t)!=0);
812  }
813 
814  if (chaine(i-1,1) == "\"") inString = !inString;
815  if (inString) continue;
816  if (chaine(i-1,1) == "[") compt2++;
817  if (chaine(i-1,1) == "]") compt2--;
818  if (chaine(i-1,1) == "(") compt++;
819  if (chaine(i-1,1) == ")") compt--;
820  if (chaine(j-1,1) == "[") compt3++;
821  if (chaine(j-1,1) == "]") compt3--;
822  if (chaine(j-1,1) == "(") compt4++;
823  if (chaine(j-1,1) == ")") compt4--;
824  if (chaine(i-1,2)=="&&" && !inString && compt==0 && compt2==0 && et==0) {et=i;puiss=0;}
825  if (chaine(i-1,2)=="||" && compt==0 && compt2==0 && ou==0) {puiss10=0; ou=i;}
826  if (chaine(i-1,1)=="&" && compt==0 && compt2==0 && etx==0) {etx=i;puiss=0;}
827  if (chaine(i-1,1)=="|" && compt==0 && compt2==0 && oux==0) {puiss10=0; oux=i;}
828  if (chaine(i-1,2)==">>" && compt==0 && compt2==0 && rshift==0) {puiss10=0; rshift=i;}
829  if (chaine(i-1,1)==">" && compt==0 && compt2==0 && rshift==0 && grand==0)
830  {puiss10=0; grand=i;}
831  if (chaine(i-1,2)=="<<" && compt==0 && compt2==0 && lshift==0) {puiss10=0; lshift=i;}
832  if (chaine(i-1,1)=="<" && compt==0 && compt2==0 && lshift==0 && petit==0)
833  {puiss10=0; petit=i;
834  // Check whether or not we have a template names! (actually this can
835  // only happen in TTreeFormula.
836  for(int ip = i,depth=0; ip < lchain; ++ip) {
837  char c = chaine(ip);
838  // The characters allowed in the template parameter are alpha-numerical characters,
839  // underscores, comma, <, > and scope operator.
840  if (isalnum(c) || c=='_' || c==',') continue;
841  if (c==':' && chaine(ip+1)==':') { ++ip; continue; }
842  if (c=='<') { ++depth; continue; }
843  if (c=='>') {
844  if (depth) { --depth; continue; }
845  else {
846  // We reach the end of the template parameter.
847  petit = 0;
848  i = ip+1;
849  break;
850  }
851  }
852  // Character not authorized within a template parameter
853  break;
854  }
855  if (petit==0) {
856  // We found a template parameter and modified i
857  continue; // the for(int i ,...)
858  }
859  }
860  if ((chaine(i-1,2)=="<=" || chaine(i-1,2)=="=<") && compt==0 && compt2==0
861  && peteg==0) {peteg=i; puiss10=0; petit=0;}
862  if ((chaine(i-1,2)=="=>" || chaine(i-1,2)==">=") && compt==0 && compt2==0
863  && grdeg==0) {puiss10=0; grdeg=i; grand=0;}
864  if (chaine(i-1,2) == "==" && compt == 0 && compt2 == 0 && egal == 0) {puiss10=0; egal=i;}
865  if (chaine(i-1,2) == "!=" && compt == 0 && compt2 == 0 && diff == 0) {puiss10=0; diff=i;}
866  if (i>1 && chaine(i-1,1) == "+" && compt == 0 && compt2 == 0 && puiss10==0) plus=i;
867  if (chaine(j-1,1) == "-" && chaine(j-2,1) != "*" && chaine(j-2,1) != "/"
868  && chaine(j-2,1)!="^" && compt3==0 && compt4==0 && moins==0 && puiss10bis==0) moins=j;
869  if (chaine(i-1,1)=="%" && compt==0 && compt2==0 && modulo==0) {puiss10=0; modulo=i;}
870  if (chaine(i-1,1)=="*" && compt==0 && compt2==0 && multi==0) {puiss10=0; multi=i;}
871  if (chaine(j-1,1)=="/" && chaine(j-2,1)!="\\"
872  && compt4==0 && compt3==0 && divi==0)
873  {
874  puiss10=0; divi=j;
875  }
876  if (chaine(j-1)=='^' && compt4==0 && compt3==0 && puiss==0) {puiss10=0; puiss=j;}
877  if (chaine(i-1)=='?' && compt == 0 && compt2 == 0 && tercond == 0) {puiss10=0; tercond=i;}
878  if (chaine(i-1)==':' && tercond && compt == 0 && compt2 == 0 && terelse == 0) {
879  if (i>2 && chaine(i-2)!=':' && chaine(i)!=':') {
880  puiss10=0; terelse=i;
881  }
882  }
883 
884  j--;
885  }
886 
887  // If operator found, analyze left and right part of the statement
888 
889  actionParam = 0;
890  if (tercond && terelse) {
891  if (tercond == 1 || terelse == lchain || tercond == (terelse-1) ) {
892  err = 5;
893  chaine_error = "?:";
894  } else {
895  // Condition
896  ctemp = chaine(0,tercond-1);
897  Analyze(ctemp.Data(),err,offset); if (err) return;
898 
899  fExpr[fNoper] = "?: condition jump";
900  actionCode = kJumpIf;
901  actionParam = 0;
902  SetAction(fNoper,actionCode, actionParam);
903  Int_t optloc = fNoper++;
904 
905  // Expression executed if condition is true.
906  ctemp = chaine(tercond,terelse-tercond-1);
907  Analyze(ctemp.Data(),err,offset); if (err) return;
908  actionParam = fNoper; // We want to skip the next instruction ('else jump'), so we set the param to the current cursor and the next instruction will be skip by the ++i in the eval loop
909  SetAction(optloc, actionCode, actionParam);
910 
911  fExpr[fNoper] = "?: else jump";
912  actionCode = kJump;
913  actionParam = 0;
914  // Set jump target.
915  SetAction(fNoper,actionCode, actionParam);
916  optloc = fNoper++;
917 
918  // Expression executed if condition is false.
919  ctemp = chaine(terelse,lchain-terelse);
920  Analyze(ctemp.Data(),err,offset); if (err) return;
921  // Set jump target.
922  actionParam = fNoper - 1; // We need to not skip the next instruction, so we compensate for the ++i in the eval loop
923  SetAction(optloc, actionCode, actionParam);
924 
925  if (IsString(optloc-1) != IsString(fNoper-1)) {
926  err = 45;
927  chaine_error = "?:";
928  }
929  }
930  } else if (ou != 0) { //check for ||
931  if (ou==1 || ou==lchain-1) {
932  err=5;
933  chaine_error="||";
934  }
935  else {
936  ctemp = chaine(0,ou-1);
937  Analyze(ctemp.Data(),err,offset); if (err) return;
938 
939  fExpr[fNoper] = "|| checkpoint";
940  actionCode = kBoolOptimize;
941  actionParam = 2;
942  SetAction(fNoper,actionCode, actionParam);
943  Int_t optloc = fNoper++;
944 
945  ctemp = chaine(ou+1,lchain-ou-1);
946  Analyze(ctemp.Data(),err,offset); if (err) return;
947  fExpr[fNoper] = "||";
948  actionCode = kOr;
949  SetAction(fNoper,actionCode, 0);
950 
951  SetAction( optloc, GetAction(optloc), GetActionParam(optloc) + (fNoper-optloc) * 10);
952  fNoper++;
953  if (!CheckOperands(optloc-1,fNoper-1,err)) return;
954  }
955  } else if (et!=0) {
956  if (et==1 || et==lchain-1) {
957  err=5;
958  chaine_error="&&";
959  }
960  else {
961  ctemp = chaine(0,et-1);
962  Analyze(ctemp.Data(),err,offset); if (err) return;
963 
964  fExpr[fNoper] = "&& checkpoint";
965  actionCode = kBoolOptimize;
966  actionParam = 1;
967  SetAction(fNoper,actionCode,actionParam);
968 
969  Int_t optloc = fNoper++;
970 
971  ctemp = chaine(et+1,lchain-et-1);
972  Analyze(ctemp.Data(),err,offset); if (err) return;
973  fExpr[fNoper] = "&&";
974  actionCode = kAnd;
975  SetAction(fNoper,actionCode,0);
976 
977  SetAction(optloc, GetAction(optloc), GetActionParam(optloc) + (fNoper-optloc) * 10);
978  fNoper++;
979  if (!CheckOperands(optloc-1,fNoper-1,err)) return;
980  }
981  } else if (oux!=0) {
982  if (oux==1 || oux==lchain) {
983  err=5;
984  chaine_error="|";
985  }
986  else {
987  ctemp = chaine(0,oux-1);
988  Analyze(ctemp.Data(),err,offset); if (err) return;
989  UInt_t leftopr = fNoper-1;
990  ctemp = chaine(oux,lchain-oux);
991  Analyze(ctemp.Data(),err,offset); if (err) return;
992  fExpr[fNoper] = "|";
993  actionCode = kBitOr;
994  SetAction(fNoper,actionCode,actionParam);
995  fNoper++;
996  if (!CheckOperands(leftopr,fNoper-1,err)) return;
997  }
998  } else if (etx!=0) {
999  if (etx==1 || etx==lchain) {
1000  err=5;
1001  chaine_error="&";
1002  }
1003  else {
1004  ctemp = chaine(0,etx-1);
1005  Analyze(ctemp.Data(),err,offset); if (err) return;
1006  UInt_t leftopr = fNoper-1;
1007  ctemp = chaine(etx,lchain-etx);
1008  Analyze(ctemp.Data(),err,offset); if (err) return;
1009  fExpr[fNoper] = "&";
1010  actionCode = kBitAnd;
1011  SetAction(fNoper,actionCode,actionParam);
1012  fNoper++;
1013  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1014  }
1015  } else if (petit != 0) {
1016  if (petit==1 || petit==lchain) {
1017  err=5;
1018  chaine_error="<";
1019  }
1020  else {
1021  ctemp = chaine(0,petit-1);
1022  Analyze(ctemp.Data(),err,offset); if (err) return;
1023  UInt_t leftopr = fNoper-1;
1024  ctemp = chaine(petit,lchain-petit);
1025  Analyze(ctemp.Data(),err,offset); if (err) return;
1026  fExpr[fNoper] = "<";
1027  actionCode = kLess;
1028  SetAction(fNoper,actionCode,actionParam);
1029  fNoper++;
1030  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1031  }
1032  } else if (grand != 0) {
1033  if (grand==1 || grand==lchain) {
1034  err=5;
1035  chaine_error=">";
1036  }
1037  else {
1038  ctemp = chaine(0,grand-1);
1039  Analyze(ctemp.Data(),err,offset); if (err) return;
1040  UInt_t leftopr = fNoper-1;
1041  ctemp = chaine(grand,lchain-grand);
1042  Analyze(ctemp.Data(),err,offset); if (err) return;
1043  fExpr[fNoper] = ">";
1044  actionCode = kGreater;
1045  SetAction(fNoper,actionCode,actionParam);
1046  fNoper++;
1047  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1048  }
1049  } else if (peteg != 0) {
1050  if (peteg==1 || peteg==lchain-1) {
1051  err=5;
1052  chaine_error="<=";
1053  }
1054  else {
1055  ctemp = chaine(0,peteg-1);
1056  Analyze(ctemp.Data(),err,offset); if (err) return;
1057  ctemp = chaine(peteg+1,lchain-peteg-1);
1058  UInt_t leftopr = fNoper-1;
1059  Analyze(ctemp.Data(),err,offset); if (err) return;
1060  fExpr[fNoper] = "<=";
1061  actionCode = kLessThan;
1062  SetAction(fNoper,actionCode,actionParam);
1063  fNoper++;
1064  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1065  }
1066  } else if (grdeg != 0) {
1067  if (grdeg==1 || grdeg==lchain-1) {
1068  err=5;
1069  chaine_error="=>";
1070  }
1071  else {
1072  ctemp = chaine(0,grdeg-1);
1073  Analyze(ctemp.Data(),err,offset); if (err) return;
1074  UInt_t leftopr = fNoper-1;
1075  ctemp = chaine(grdeg+1,lchain-grdeg-1);
1076  Analyze(ctemp.Data(),err,offset); if (err) return;
1077  fExpr[fNoper] = ">=";
1078  actionCode = kGreaterThan;
1079  SetAction(fNoper,actionCode,actionParam);
1080  fNoper++;
1081  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1082  }
1083  } else if (egal != 0) {
1084  if (egal==1 || egal==lchain-1) {
1085  err=5;
1086  chaine_error="==";
1087  }
1088  else {
1089  ctemp = chaine(0,egal-1);
1090  Analyze(ctemp.Data(),err,offset); if (err) return;
1091  Int_t optloc = fNoper-1;
1092 
1093  ctemp = chaine(egal+1,lchain-egal-1);
1094  Analyze(ctemp.Data(),err,offset); if (err) return;
1095  fExpr[fNoper] = "==";
1096  actionCode = kEqual;
1097 
1098  Bool_t isstring = IsString(fNoper-1);
1099  if (IsString(optloc) != isstring) {
1100  err = 45;
1101  chaine_error = "==";
1102  } else if (isstring) {
1103  actionCode = kStringEqual;
1104  }
1105  SetAction(fNoper,actionCode,actionParam);
1106  fNoper++;
1107  }
1108  } else if (diff != 0) {
1109  if (diff==1 || diff==lchain-1) {
1110  err=5;
1111  chaine_error = "!=";
1112  }
1113  else {
1114  ctemp = chaine(0,diff-1);
1115  Analyze(ctemp.Data(),err,offset); if (err) return;
1116  Int_t optloc = fNoper-1;
1117 
1118  ctemp = chaine(diff+1,lchain-diff-1);
1119  Analyze(ctemp.Data(),err,offset); if (err) return;
1120  fExpr[fNoper] = "!=";
1121  actionCode = kNotEqual;
1122 
1123  Bool_t isstring = IsString(fNoper-1);
1124  if (IsString(optloc) != isstring) {
1125  err = 45;
1126  chaine_error = "!=";
1127  } else if (isstring) {
1128  actionCode = kStringNotEqual;
1129  }
1130  SetAction(fNoper,actionCode,actionParam);
1131  fNoper++;
1132  }
1133  } else if (plus != 0) {
1134  if (plus==lchain) {
1135  err=5;
1136  chaine_error = "+";
1137  }
1138  else {
1139  ctemp = chaine(0,plus-1);
1140  Analyze(ctemp.Data(),err,offset); if (err) return;
1141  UInt_t leftopr = fNoper-1;
1142  ctemp = chaine(plus,lchain-plus);
1143  Analyze(ctemp.Data(),err,offset); if (err) return;
1144  fExpr[fNoper] = "+";
1145  actionCode = kAdd;
1146  SetAction(fNoper,actionCode,actionParam);
1147  fNoper++;
1148  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1149  }
1150  } else {
1151  if (moins != 0) {
1152  if (moins == 1) {
1153  ctemp = chaine(moins,lchain-moins);
1154  Analyze(ctemp.Data(),err,offset); if (err) return;
1155  fExpr[fNoper] = "-";
1156  actionCode = kSignInv;
1157  SetAction(fNoper,actionCode,actionParam);
1158  ++fNoper;
1159  if (!CheckOperands(fNoper-1,err)) return;
1160  } else {
1161  if (moins == lchain) {
1162  err=5;
1163  chaine_error = "-";
1164  } else {
1165  ctemp = chaine(0,moins-1);
1166  Analyze(ctemp.Data(),err,offset); if (err) return;
1167  UInt_t leftopr = fNoper-1;
1168  ctemp = chaine(moins,lchain-moins);
1169  Analyze(ctemp.Data(),err,offset); if (err) return;
1170  fExpr[fNoper] = "-";
1171  actionCode = kSubstract;
1172  SetAction(fNoper,actionCode,actionParam);
1173  fNoper++;
1174  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1175  }
1176  }
1177  } else if (modulo != 0) {
1178  if (modulo == 1 || modulo == lchain) {
1179  err=5;
1180  chaine_error="%";
1181  } else {
1182  ctemp = chaine(0,modulo-1);
1183  Analyze(ctemp.Data(),err,offset); if (err) return;
1184  UInt_t leftopr = fNoper-1;
1185  ctemp = chaine(modulo,lchain-modulo);
1186  Analyze(ctemp.Data(),err,offset); if (err) return;
1187  fExpr[fNoper] = "%";
1188  actionCode = kModulo;
1189  SetAction(fNoper,actionCode,actionParam);
1190  fNoper++;
1191  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1192  }
1193  } else if (rshift != 0) {
1194  if (rshift == 1 || rshift == lchain) {
1195  err=5;
1196  chaine_error=">>";
1197  } else {
1198  ctemp = chaine(0,rshift-1);
1199  Analyze(ctemp.Data(),err,offset); if (err) return;
1200  UInt_t leftopr = fNoper-1;
1201  ctemp = chaine(rshift+1,lchain-rshift-1);
1202  Analyze(ctemp.Data(),err,offset); if (err) return;
1203  fExpr[fNoper] = ">>";
1204  actionCode = kRightShift;
1205  SetAction(fNoper,actionCode,actionParam);
1206  fNoper++;
1207  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1208  }
1209  } else if (lshift != 0) {
1210  if (lshift == 1 || lshift == lchain) {
1211  err=5;
1212  chaine_error=">>";
1213  } else {
1214  ctemp = chaine(0,lshift-1);
1215  Analyze(ctemp.Data(),err,offset); if (err) return;
1216  UInt_t leftopr = fNoper-1;
1217  ctemp = chaine(lshift+1,lchain-lshift-1);
1218  Analyze(ctemp.Data(),err,offset); if (err) return;
1219  fExpr[fNoper] = ">>";
1220  actionCode = kLeftShift;
1221  SetAction(fNoper,actionCode,actionParam);
1222  fNoper++;
1223  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1224  }
1225  } else {
1226  if (multi != 0) {
1227  if (multi == 1 || multi == lchain) {
1228  err=5;
1229  chaine_error="*";
1230  }
1231  else {
1232  ctemp = chaine(0,multi-1);
1233  Analyze(ctemp.Data(),err,offset); if (err) return;
1234  UInt_t leftopr = fNoper-1;
1235  ctemp = chaine(multi,lchain-multi);
1236  Analyze(ctemp.Data(),err,offset); if (err) return;
1237  fExpr[fNoper] = "*";
1238  actionCode = kMultiply;
1239  SetAction(fNoper,actionCode,actionParam);
1240  fNoper++;
1241  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1242  }
1243  } else {
1244  if (divi != 0) {
1245  if (divi == 1 || divi == lchain) {
1246  err=5;
1247  chaine_error = "/";
1248  }
1249  else {
1250  ctemp = chaine(0,divi-1);
1251  Analyze(ctemp.Data(),err,offset); if (err) return;
1252  UInt_t leftopr = fNoper-1;
1253  ctemp = chaine(divi,lchain-divi);
1254  Analyze(ctemp.Data(),err,offset); if (err) return;
1255  fExpr[fNoper] = "/";
1256  actionCode = kDivide;
1257  SetAction(fNoper,actionCode,actionParam);
1258  fNoper++;
1259  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1260  }
1261  } else {
1262  if (puiss != 0) {
1263  if (puiss == 1 || puiss == lchain) {
1264  err = 5;
1265  chaine_error = "**";
1266  }
1267  else {
1268  if (chaine(lchain-2,2) == "^2") {
1269  ctemp = "sq(" + chaine(0,lchain-2) + ")";
1270  Analyze(ctemp.Data(),err,offset); if (err) return;
1271  } else {
1272  ctemp = chaine(0,puiss-1);
1273  Analyze(ctemp.Data(),err,offset); if (err) return;
1274  UInt_t leftopr = fNoper-1;
1275  ctemp = chaine(puiss,lchain-puiss);
1276  Analyze(ctemp.Data(),err,offset); if (err) return;
1277  fExpr[fNoper] = "^";
1278  actionCode = kpow;
1279  SetAction(fNoper,actionCode,actionParam);
1280  fNoper++;
1281  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1282  }
1283  }
1284  } else {
1285 
1286  find=0;
1287 
1288  // Check for a numerical expression
1289  {
1290  Bool_t hasDot = kFALSE;
1291  Bool_t isHexa = kFALSE;
1292  Bool_t hasExpo= kFALSE;
1293  if ((chaine(0,2)=="0x")||(chaine(0,2)=="0X")) isHexa=kTRUE;
1294  for (j=0; j<chaine.Length() && err==0; j++) {
1295  t=chaine[j];
1296  if (!isHexa) {
1297  if (j>0 && (chaine(j,1)=="e" || chaine(j,2)=="e+" || chaine(j,2)=="e-" || chaine(j,1)=="E" || chaine(j,2)=="E+" || chaine(j,2)=="E-")) {
1298  if (hasExpo) {
1299  err=26;
1300  chaine_error=chaine;
1301  }
1302  hasExpo = kTRUE;
1303  // The previous implementation allowed a '.' in the exponent.
1304  // That information was ignored (by sscanf), we now make it an error
1305  // hasDot = kFALSE;
1306  hasDot = kTRUE; // forbid any additional '.'
1307  if (chaine(j,2)=="e+" || chaine(j,2)=="e-" || chaine(j,2)=="E+" || chaine(j,2)=="E-") j++;
1308  }
1309  else {
1310  if (chaine(j,1) == "." && !hasDot) hasDot = kTRUE; // accept only one '.' in the number
1311  else {
1312  // The previous implementation was allowing ANYTHING after the '.' and thus
1313  // made legal code like '2.3 and fpx' and was just silently ignoring the
1314  // 'and fpx'.
1315  if (!strchr("0123456789",t) && (chaine(j,1)!="+" || j!=0)) {
1316  err = 30;
1317  chaine_error=chaine;
1318  }
1319  }
1320  }
1321  }
1322  else {
1323  if (!strchr("0123456789abcdefABCDEF",t) && (j>1)) {
1324  err = 30;
1325  chaine_error=chaine;
1326  }
1327  }
1328  }
1329  if (fNconst >= gMAXCONST) err = 27;
1330  if (!err) {
1331  if (!isHexa) {if (sscanf((const char*)chaine,"%lg",&vafConst) > 0) err = 0; else err =1;}
1332  else {if (sscanf((const char*)chaine,"%lx",&vafConst2) > 0) err = 0; else err=1;
1333  vafConst = (Double_t) vafConst2;}
1334  fExpr[fNoper] = chaine;
1335  k = -1;
1336  for (j=0;j<fNconst;j++) {
1337  if (vafConst == fConst[j] ) k= j;
1338  }
1339  if ( k < 0) { k = fNconst; fNconst++; fConst[k] = vafConst; }
1340  actionCode = kConstant;
1341  actionParam = k;
1342  SetAction(fNoper,actionCode,actionParam);
1343  fNoper++;
1344  }
1345  if (err==30) err=0;
1346  else find = kTRUE;
1347  }
1348 
1349  // Look for an already defined expression
1350 
1351  if (find==0) {
1352  {
1353  R__LOCKGUARD(gROOTMutex);
1354  oldformula = (const TFormula*)gROOT->GetListOfFunctions()->FindObject((const char*)chaine);
1355  }
1356  if (oldformula && strcmp(schain,oldformula->GetTitle())) {
1357  Int_t nprior = fNpar;
1358  Analyze(oldformula->GetExpFormula(),err,fNpar);
1359  // if the oldformula was using a normalized function (gausn or landaun) set also in this one
1360  if (oldformula->IsNormalized()) SetBit(kNormalized);
1361  if (err) return; // changes fNpar
1362  fNpar = nprior;
1363  find=1;
1364  if (!err) {
1365  Int_t npold = oldformula->GetNpar();
1366  fNpar += npold;
1367  for (Int_t ipar=0;ipar<npold;ipar++) {
1368  fParams[ipar+fNpar-npold] = oldformula->GetParameter(ipar);
1369  }
1370  }
1371  }
1372  }
1373  if (find == 0) {
1374 
1375  // Check if chaine is a defined variable.
1376  // Note that DefinedVariable can be overloaded
1377 
1378  ctemp = chaine;
1379  ctemp.ReplaceAll(escapedSlash, slash);
1380  Int_t action;
1381  k = DefinedVariable(ctemp,action);
1382  if (k==-3) {
1383  // Error message already issued
1384  err = 1;
1385  } else if (k==-2) {
1386  err = 31;
1387  chaine_error = ctemp;
1388  } else if ( k >= 0 ) {
1389  fExpr[fNoper] = ctemp;
1390  actionCode = action;
1391  actionParam = k;
1392  SetAction(fNoper,actionCode,actionParam);
1393  if (action==kDefinedString) fNstring++;
1394  else if (k <kMAXFOUND && !fAlreadyFound.TestBitNumber(k)) {
1395  fAlreadyFound.SetBitNumber(k);
1396  fNval++;
1397  }
1398  fNoper++;
1399  } else if (chaine(0,1) == "!") {
1400  ctemp = chaine(1,lchain-1);
1401  Analyze(ctemp.Data(),err,offset); if (err) return;
1402  fExpr[fNoper] = "!";
1403  actionCode = kNot;
1404  SetAction(fNoper,actionCode,actionParam);
1405  fNoper++;
1406  if (!CheckOperands(fNoper-1,err)) return;
1407  } else if (chaine(0,1)=="\"" && chaine(chaine.Length()-1,1)=="\"") {
1408  // It is a string !!!
1409  fExpr[fNoper] = chaine(1,chaine.Length()-2);
1410  actionCode = kStringConst;
1411  SetAction(fNoper,actionCode,actionParam);
1412  fNoper++;
1413  } else if (chaine(0,4) == "cos(") {
1414  ctemp = chaine(3,lchain-3);
1415  Analyze(ctemp.Data(),err,offset); if (err) return;
1416  fExpr[fNoper] = "cos";
1417  actionCode = kcos;
1418  SetAction(fNoper,actionCode,actionParam);
1419  fNoper++;
1420  if (!CheckOperands(fNoper-1,err)) return;
1421  } else if (chaine(0,4) == "sin(") {
1422  ctemp = chaine(3,lchain-3);
1423  Analyze(ctemp.Data(),err,offset); if (err) return;
1424  fExpr[fNoper] = "sin";
1425  actionCode = ksin;
1426  SetAction(fNoper,actionCode,actionParam);
1427  fNoper++;
1428  if (!CheckOperands(fNoper-1,err)) return;
1429  } else if (chaine(0,4) == "tan(") {
1430  ctemp = chaine(3,lchain-3);
1431  Analyze(ctemp.Data(),err,offset); if (err) return;
1432  fExpr[fNoper] = "tan";
1433  actionCode = ktan;
1434  SetAction(fNoper,actionCode,actionParam);
1435  fNoper++;
1436  if (!CheckOperands(fNoper-1,err)) return;
1437  } else if (chaine(0,5) == "acos(") {
1438  ctemp = chaine(4,lchain-4);
1439  Analyze(ctemp.Data(),err,offset); if (err) return;
1440  fExpr[fNoper] = "acos";
1441  actionCode = kacos;
1442  SetAction(fNoper,actionCode,actionParam);
1443  fNoper++;
1444  if (!CheckOperands(fNoper-1,err)) return;
1445  } else if (chaine(0,5) == "asin(") {
1446  ctemp = chaine(4,lchain-4);
1447  Analyze(ctemp.Data(),err,offset); if (err) return;
1448  fExpr[fNoper] = "asin";
1449  actionCode = kasin;
1450  SetAction(fNoper,actionCode,actionParam);
1451  fNoper++;
1452  if (!CheckOperands(fNoper-1,err)) return;
1453  } else if (chaine(0,5) == "atan(") {
1454  ctemp = chaine(4,lchain-4);
1455  Analyze(ctemp.Data(),err,offset); if (err) return;
1456  fExpr[fNoper] = "atan";
1457  actionCode = katan;
1458  SetAction(fNoper,actionCode,actionParam);
1459  fNoper++;
1460  if (!CheckOperands(fNoper-1,err)) return;
1461  } else if (chaine(0,5) == "cosh(") {
1462  ctemp = chaine(4,lchain-4);
1463  Analyze(ctemp.Data(),err,offset); if (err) return;
1464  fExpr[fNoper] = "cosh";
1465  actionCode = kcosh;
1466  SetAction(fNoper,actionCode,actionParam);
1467  fNoper++;
1468  if (!CheckOperands(fNoper-1,err)) return;
1469  } else if (chaine(0,5) == "sinh(") {
1470  ctemp = chaine(4,lchain-4);
1471  Analyze(ctemp.Data(),err,offset); if (err) return;
1472  fExpr[fNoper] = "sinh";
1473  actionCode = ksinh;
1474  SetAction(fNoper,actionCode,actionParam);
1475  fNoper++;;
1476  if (!CheckOperands(fNoper-1,err)) return;
1477  } else if (chaine(0,5) == "tanh(") {
1478  ctemp = chaine(4,lchain-4);
1479  Analyze(ctemp.Data(),err,offset); if (err) return;
1480  fExpr[fNoper] = "tanh";
1481  actionCode = ktanh;
1482  SetAction(fNoper,actionCode,actionParam);
1483  fNoper++;;
1484  if (!CheckOperands(fNoper-1,err)) return;
1485  } else if (chaine(0,6) == "acosh(") {
1486  ctemp = chaine(5,lchain-5);
1487  Analyze(ctemp.Data(),err,offset); if (err) return;
1488  fExpr[fNoper] = "acosh";
1489  actionCode = kacosh;
1490  SetAction(fNoper,actionCode,actionParam);
1491  fNoper++;;
1492  if (!CheckOperands(fNoper-1,err)) return;
1493  } else if (chaine(0,6) == "asinh(") {
1494  ctemp = chaine(5,lchain-5);
1495  Analyze(ctemp.Data(),err,offset); if (err) return;
1496  fExpr[fNoper] = "asinh";
1497  actionCode = kasinh;
1498  SetAction(fNoper,actionCode,actionParam);
1499  fNoper++;;
1500  if (!CheckOperands(fNoper-1,err)) return;
1501  } else if (chaine(0,6) == "atanh(") {
1502  ctemp = chaine(5,lchain-5);
1503  Analyze(ctemp.Data(),err,offset); if (err) return;
1504  fExpr[fNoper] = "atanh";
1505  actionCode = katanh;
1506  SetAction(fNoper,actionCode,actionParam);
1507  fNoper++;;
1508  if (!CheckOperands(fNoper-1,err)) return;
1509  } else if (chaine(0,3) == "sq(") {
1510  ctemp = chaine(2,lchain-2);
1511  Analyze(ctemp.Data(),err,offset); if (err) return;
1512  fExpr[fNoper] = "sq";
1513  actionCode = ksq;
1514  SetAction(fNoper,actionCode,actionParam);
1515  fNoper++;;
1516  if (!CheckOperands(fNoper-1,err)) return;
1517  } else if (chaine(0,4) == "log(") {
1518  ctemp = chaine(3,lchain-3);
1519  Analyze(ctemp.Data(),err,offset); if (err) return;
1520  fExpr[fNoper] = "log";
1521  actionCode = klog;
1522  SetAction(fNoper,actionCode,actionParam);
1523  fNoper++;;
1524  if (!CheckOperands(fNoper-1,err)) return;
1525  } else if (chaine(0,6) == "log10(") {
1526  ctemp = chaine(5,lchain-5);
1527  Analyze(ctemp.Data(),err,offset); if (err) return;
1528  fExpr[fNoper] = "log10";
1529  actionCode = klog10;
1530  SetAction(fNoper,actionCode,actionParam);
1531  fNoper++;;
1532  if (!CheckOperands(fNoper-1,err)) return;
1533  } else if (chaine(0,4) == "exp(") {
1534  ctemp = chaine(3,lchain-3);
1535  Analyze(ctemp.Data(),err,offset); if (err) return;
1536  fExpr[fNoper] = "exp";
1537  actionCode = kexp;
1538  SetAction(fNoper,actionCode,actionParam);
1539  fNoper++;;
1540  if (!CheckOperands(fNoper-1,err)) return;
1541  } else if (chaine(0,4) == "abs(") {
1542  ctemp = chaine(3,lchain-3);
1543  Analyze(ctemp.Data(),err,offset); if (err) return;
1544  fExpr[fNoper] = "abs";
1545  actionCode = kabs;
1546  SetAction(fNoper,actionCode,actionParam);
1547  fNoper++;;
1548  if (!CheckOperands(fNoper-1,err)) return;
1549  } else if (chaine(0,5) == "sign(") {
1550  ctemp = chaine(4,lchain-4);
1551  Analyze(ctemp.Data(),err,offset); if (err) return;
1552  fExpr[fNoper] = "sign";
1553  actionCode = ksign;
1554  SetAction(fNoper,actionCode,actionParam);
1555  fNoper++;;
1556  if (!CheckOperands(fNoper-1,err)) return;
1557  } else if (chaine(0,4) == "int(") {
1558  ctemp = chaine(3,lchain-3);
1559  Analyze(ctemp.Data(),err,offset); if (err) return;
1560  fExpr[fNoper] = "int";
1561  actionCode = kint;
1562  SetAction(fNoper,actionCode,actionParam);
1563  fNoper++;;
1564  if (!CheckOperands(fNoper-1,err)) return;
1565  } else if (chaine == "rndm" || chaine(0,5) == "rndm(") {
1566  fExpr[fNoper] = "rndm";
1567  actionCode = krndm;
1568  SetAction(fNoper,actionCode,actionParam);
1569  fNoper++;;
1570  } else if (chaine(0,5) == "sqrt(") {
1571  ctemp = chaine(4,lchain-4);
1572  Analyze(ctemp.Data(),err,offset); if (err) return;
1573  fExpr[fNoper] = "sqrt";
1574  actionCode = ksqrt;
1575  SetAction(fNoper,actionCode,actionParam);
1576  fNoper++;;
1577  if (!CheckOperands(fNoper-1,err)) return;
1578 
1579  // Look for an exponential
1580 
1581  } else if ( chaine == "expo" || chaine(0,5)=="expo("
1582  || (lchain==5 && chaine(1,4)=="expo")
1583  || (lchain==6 && chaine(2,4)=="expo")
1584  || chaine(1,5)=="expo(" || chaine(2,5)=="expo(" ) {
1585  chaine1ST=chaine;
1586  if (chaine(1,4) == "expo") {
1587  ctemp=chaine(0,1);
1588  if (ctemp=="x") {
1589  inter2=0;
1590  if (fNdim < 1) fNdim = 1; }
1591  else if (ctemp=="y") {
1592  inter2=1;
1593  if (fNdim < 2) fNdim = 2; }
1594  else if (ctemp=="z") {
1595  inter2=2;
1596  if (fNdim < 3) fNdim = 3; }
1597  else if (ctemp=="t") {
1598  inter2=3;
1599  if (fNdim < 4) fNdim = 4; }
1600  else {
1601  err=26; // unknown name;
1602  chaine_error=chaine1ST;
1603  }
1604  chaine=chaine(1,lchain-1);
1605  lchain=chaine.Length();
1606  } else inter2=0;
1607  if (chaine(2,4) == "expo") {
1608  if (chaine(0,2) != "xy") {
1609  err=26; // unknown name
1610  chaine_error=chaine1ST;
1611  }
1612  else {
1613  inter2=5;
1614  if (fNdim < 2) fNdim = 2;
1615  chaine=chaine(2,lchain-2);
1616  lchain=chaine.Length();
1617  }
1618  }
1619  if (lchain == 4) {
1620  if (fNpar>=gMAXPAR) err=7; // too many parameters
1621  if (!err) {
1622  fExpr[fNoper] = chaine1ST;
1623  actionCode = kexpo + inter2;
1624  actionParam = offset;
1625  SetAction(fNoper,actionCode,actionParam);
1626  if (inter2 == 5+offset && fNpar < 3+offset) fNpar = 3+offset;
1627  if (fNpar < 2+offset) fNpar = 2+offset;
1628  if (fNpar>=gMAXPAR) err=7; // too many parameters
1629  if (!err) {
1630  fNoper++;
1631  if (fNdim < 1) fNdim = 1;
1632  if (fNpar == 2) SetNumber(200);
1633  }
1634  }
1635  } else if (chaine(4,1) == "(") {
1636  ctemp = chaine(5,lchain-6);
1637  fExpr[fNoper] = chaine1ST;
1638  for (j=0; j<ctemp.Length(); j++) {
1639  t=ctemp[j];
1640  if (strchr("0123456789",t)==0 && (ctemp(j,1)!="+" || j!=0)) {
1641  err=20;
1642  chaine_error=chaine1ST;
1643  }
1644  }
1645  if (err==0) {
1646  sscanf(ctemp.Data(),"%d",&inter);
1647  if (inter>=0) {
1648  inter += offset;
1649  actionCode = kexpo + inter2;
1650  actionParam = inter;
1651  SetAction(fNoper,actionCode,actionParam);
1652  if (inter2 == 5) inter++;
1653  if (inter+2>fNpar) fNpar = inter+2;
1654  if (fNpar>=gMAXPAR) err=7; // too many parameters
1655  if (!err) fNoper++;
1656  if (fNpar == 2) SetNumber(200);
1657  } else err=20;
1658  } else err = 20; // non integer value for parameter number
1659  } else {
1660  err=26; // unknown name
1661  chaine_error=chaine;
1662  }
1663 
1664  // Look for gaus, xgaus,ygaus,xygaus
1665 
1666  } else if (chaine=="gaus"
1667  || (lchain==5 && chaine(1,4)=="gaus")
1668  || (lchain==6 && chaine(2,4)=="gaus")
1669  || chaine(0,5)=="gaus(" || chaine(1,5)=="gaus(" || chaine(2,5)=="gaus(") {
1670  chaine1ST=chaine;
1671  if (chaine(1,4) == "gaus") {
1672  ctemp=chaine(0,1);
1673  if (ctemp=="x") {
1674  inter2=0;
1675  if (fNdim < 1) fNdim = 1; }
1676  else if (ctemp=="y") {
1677  inter2=1;
1678  if (fNdim < 2) fNdim = 2; }
1679  else if (ctemp=="z") {
1680  inter2=2;
1681  if (fNdim < 3) fNdim = 3; }
1682  else if (ctemp=="t") {
1683  inter2=3;
1684  if (fNdim < 4) fNdim = 4; }
1685  else {
1686  err=26; // unknown name
1687  chaine_error=chaine1ST;
1688  }
1689  chaine=chaine(1,lchain-1);
1690  lchain=chaine.Length();
1691  } else inter2=0;
1692  if (chaine(2,4) == "gaus") {
1693  if (chaine(0,2) != "xy") {
1694  err=26; // unknown name
1695  chaine_error=chaine1ST;
1696  }
1697  else {
1698  inter2=5;
1699  if (fNdim < 2) fNdim = 2;
1700  chaine=chaine(2,lchain-2);
1701  lchain=chaine.Length();
1702  SetNumber(110); // xygaus
1703  }
1704  }
1705  if (lchain == 4 && err==0) {
1706  if (fNpar>=gMAXPAR) err=7; // too many parameters
1707  if (!err) {
1708  fExpr[fNoper] = chaine1ST;
1709  actionCode = kgaus + inter2;
1710  actionParam = offset;
1711  SetAction(fNoper,actionCode,actionParam);
1712  if (inter2 == 5+offset && fNpar < 5+offset) fNpar = 5+offset;
1713  if (3+offset>fNpar) fNpar = 3+offset;
1714  if (fNpar>=gMAXPAR) err=7; // too many parameters
1715  if (!err) {
1716  fNoper++;
1717  if (fNdim < 1) fNdim = 1;
1718  if (fNpar == 3) SetNumber(100);
1719  }
1720  }
1721  } else if (chaine(4,1) == "(" && err==0) {
1722  ctemp = chaine(5,lchain-6);
1723  fExpr[fNoper] = chaine1ST;
1724  for (j=0; j<ctemp.Length(); j++) {
1725  t=ctemp[j];
1726  if (strchr("0123456789",t)==0 && (ctemp(j,1)!="+" || j!=0)) {
1727  err=20;
1728  chaine_error=chaine1ST;
1729  }
1730  }
1731  if (err==0) {
1732  sscanf(ctemp.Data(),"%d",&inter);
1733  if (inter >= 0) {
1734  inter += offset;
1735  actionCode = kgaus + inter2;
1736  actionParam = inter;
1737  SetAction(fNoper,actionCode,actionParam);
1738  if (inter2 == 5) inter += 2;
1739  if (inter+3>fNpar) fNpar = inter+3;
1740  if (fNpar>=gMAXPAR) err=7; // too many parameters
1741  if (!err) fNoper++;
1742  if(fNpar == 3) SetNumber(100);
1743  } else err = 20; // non integer value for parameter number
1744  }
1745  } else if (err==0) {
1746  err=26; // unknown name
1747  chaine_error=chaine1ST;
1748  }
1749 
1750  // Look for landau, xlandau,ylandau,xylandau
1751 
1752  } else if (chaine=="landau" || (lchain==7 && chaine(1,6)=="landau")
1753  || (lchain==8 && chaine(2,6)=="landau")
1754  || chaine(0,7)=="landau(" || chaine(1,7)=="landau(" || chaine(2,7)=="landau(") {
1755  chaine1ST=chaine;
1756  if (chaine(1,6) == "landau") {
1757  ctemp=chaine(0,1);
1758  if (ctemp=="x") {
1759  inter2=0;
1760  if (fNdim < 1) fNdim = 1; }
1761  else if (ctemp=="y") {
1762  inter2=1;
1763  if (fNdim < 2) fNdim = 2; }
1764  else if (ctemp=="z") {
1765  inter2=2;
1766  if (fNdim < 3) fNdim = 3; }
1767  else if (ctemp=="t") {
1768  inter2=3;
1769  if (fNdim < 4) fNdim = 4; }
1770  else {
1771  err=26; // unknown name
1772  chaine_error=chaine1ST;
1773  }
1774  chaine=chaine(1,lchain-1);
1775  lchain=chaine.Length();
1776  } else inter2=0;
1777  if (chaine(2,6) == "landau") {
1778  if (chaine(0,2) != "xy") {
1779  err=26; // unknown name
1780  chaine_error=chaine1ST;
1781  }
1782  else {
1783  inter2=5;
1784  if (fNdim < 2) fNdim = 2;
1785  chaine=chaine(2,lchain-2);
1786  lchain=chaine.Length();
1787  SetNumber(410);
1788  }
1789  }
1790  if (lchain == 6 && err==0) {
1791  if (fNpar>=gMAXPAR) err=7; // too many parameters
1792  if (!err) {
1793  fExpr[fNoper] = chaine1ST;
1794  actionCode = klandau + inter2;
1795  actionParam = offset;
1796  SetAction(fNoper,actionCode,actionParam);
1797  if (inter2 == 5+offset && fNpar < 5+offset) fNpar = 5+offset;
1798  if (3+offset>fNpar) fNpar = 3+offset;
1799  if (fNpar>=gMAXPAR) err=7; // too many parameters
1800  if (!err) {
1801  fNoper++;
1802  if (fNdim < 1) fNdim = 1;
1803  if (fNpar == 3) SetNumber(400);
1804  }
1805  }
1806  } else if (chaine(6,1) == "(" && err==0) {
1807  ctemp = chaine(7,lchain-8);
1808  fExpr[fNoper] = chaine1ST;
1809  for (j=0; j<ctemp.Length(); j++) {
1810  t=ctemp[j];
1811  if (strchr("0123456789",t)==0 && (ctemp(j,1)!="+" || j!=0)) {
1812  err=20;
1813  chaine_error=chaine1ST;
1814  }
1815  }
1816  if (err==0) {
1817  sscanf(ctemp.Data(),"%d",&inter);
1818  if (inter >= 0) {
1819  inter += offset;
1820  actionCode = klandau + inter2;
1821  actionParam = inter;
1822  SetAction(fNoper,actionCode,actionParam);
1823  if (inter2 == 5) inter += 2;
1824  if (inter+3>fNpar) fNpar = inter+3;
1825  if (fNpar>=gMAXPAR) err=7; // too many parameters
1826  if (!err) fNoper++;
1827  if (fNpar == 3) SetNumber(400);
1828  } else err = 20; // non integer value for parameter number
1829  }
1830  } else if (err==0) {
1831  err=26; // unknown name
1832  chaine_error=chaine1ST;
1833  }
1834 
1835  // Look for a polynomial
1836 
1837  } else if (chaine(0,3) == "pol" || chaine(1,3) == "pol") {
1838  chaine1ST=chaine;
1839  if (chaine(1,3) == "pol") {
1840  ctemp=chaine(0,1);
1841  if (ctemp=="x") {
1842  inter2=1;
1843  if (fNdim < 1) fNdim = 1; }
1844  else if (ctemp=="y") {
1845  inter2=2;
1846  if (fNdim < 2) fNdim = 2; }
1847  else if (ctemp=="z") {
1848  inter2=3;
1849  if (fNdim < 3) fNdim = 3; }
1850  else if (ctemp=="t") {
1851  inter2=4;
1852  if (fNdim < 4) fNdim = 4; }
1853  else {
1854  err=26; // unknown name;
1855  chaine_error=chaine1ST;
1856  }
1857  chaine=chaine(1,lchain-1);
1858  lchain=chaine.Length();
1859  } else inter2=1;
1860  if (chaine(lchain-1,1) == ")") {
1861  nomb = 0;
1862  for (j=3;j<lchain;j++) if (chaine(j,1)=="(" && nomb == 0) nomb = j;
1863  if (nomb == 3) err = 23; // degree of polynomial not specified
1864  if (nomb == 0) err = 40; // '(' is expected
1865  ctemp = chaine(nomb+1,lchain-nomb-2);
1866  for (j=0; j<ctemp.Length(); j++) {
1867  t=ctemp[j];
1868  if (strchr("0123456789",t)==0 && (ctemp(j,1)!="+" || j!=0)) {
1869  err=20;
1870  chaine_error=chaine1ST;
1871  }
1872  }
1873  if (!err) {
1874  sscanf(ctemp.Data(),"%d",&inter);
1875  if (inter < 0) err = 20;
1876  }
1877  }
1878  else {
1879  nomb = lchain;
1880  inter = 0;
1881  }
1882  if (!err) {
1883  inter--;
1884  ctemp = chaine(3,nomb-3);
1885  if (sscanf(ctemp.Data(),"%d",&n) > 0) {
1886  if (n < 0 ) err = 24; //Degree of polynomial must be positive
1887  if (n >= 20) err = 25; //Degree of polynomial must be less than 20
1888  } else err = 20;
1889  }
1890  if (!err) {
1891  fExpr[fNoper] = chaine1ST;
1892  actionCode = kpol+(inter2-1);
1893  actionParam = n*100+inter+2;
1894  SetAction(fNoper,actionCode,actionParam);
1895  if (inter+n+1>=fNpar) fNpar = inter + n + 2;
1896  if (fNpar>=gMAXPAR) err=7; // too many parameters
1897  if (!err) {
1898  fNoper++;
1899  if (fNdim < 1) fNdim = 1;
1900  SetNumber(300+n);
1901  }
1902  }
1903 
1904  // Look for pow,atan2,etc
1905 
1906  } else if (chaine(0,4) == "pow(") {
1907  compt = 4; nomb = 0; virgule = 0; nest=0;
1908  while(compt != lchain) {
1909  compt++;
1910  if (chaine(compt-1,1) == "(") nest++;
1911  else if (chaine(compt-1,1) == ")") nest--;
1912  else if (chaine(compt-1,1) == "," && nest==0) {
1913  nomb++;
1914  if (nomb == 1 && virgule == 0) virgule = compt;
1915  }
1916  }
1917  if (nomb != 1) err = 22; // There are plus or minus than 2 arguments for pow
1918  else {
1919  ctemp = chaine(4,virgule-5);
1920  Analyze(ctemp.Data(),err,offset); if (err) return;
1921  UInt_t leftopr = fNoper-1;
1922  ctemp = chaine(virgule,lchain-virgule-1);
1923  Analyze(ctemp.Data(),err,offset); if (err) return;
1924  fExpr[fNoper] = "^";
1925  actionCode = kpow;
1926  SetAction(fNoper,actionCode,actionParam);
1927  fNoper++;
1928  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1929  }
1930  } else if (chaine(0,7) == "strstr(") {
1931  compt = 7; nomb = 0; virgule = 0; nest=0;
1932  inString = false;
1933  while(compt != lchain) {
1934  compt++;
1935  if (chaine(compt-1,1) == "\"") {
1936  inString = !inString;
1937  } else if (!inString) {
1938  if (chaine(compt-1,1) == "(") nest++;
1939  else if (chaine(compt-1,1) == ")") nest--;
1940  else if (chaine(compt-1,1) == "," && nest==0) {
1941  nomb++;
1942  if (nomb == 1 && virgule == 0) virgule = compt;
1943  }
1944  }
1945  }
1946  if (nomb != 1) err = 28; // There are plus or minus than 2 arguments for strstr
1947  else {
1948  ctemp = chaine(7,virgule-8);
1949  Analyze(ctemp.Data(),err,offset); if (err) return;
1950  Int_t optloc = fNoper-1;
1951 
1952  ctemp = chaine(virgule,lchain-virgule-1);
1953  Analyze(ctemp.Data(),err,offset); if (err) return;
1954  fExpr[fNoper] = "strstr";
1955  actionCode = kstrstr;
1956  SetAction(fNoper,actionCode,actionParam);
1957  fNoper++;
1958 
1959  if ( !IsString(optloc) || !IsString(fNoper-2) ) {
1960  err = 46;
1961  chaine_error = "strstr";
1962  }
1963  }
1964  } else if (chaine(0,4) == "min(") {
1965  compt = 4; nomb = 0; virgule = 0; nest=0;
1966  while(compt != lchain) {
1967  compt++;
1968  if (chaine(compt-1,1) == "(") nest++;
1969  else if (chaine(compt-1,1) == ")") nest--;
1970  else if (chaine(compt-1,1) == "," && nest==0) {
1971  nomb++;
1972  if (nomb == 1 && virgule == 0) virgule = compt;
1973  }
1974  }
1975  if (nomb != 1) {
1976  err = 44; // There are plus or minus than 2 arguments for min
1977  err_hint = 3;
1978  }
1979  else {
1980  ctemp = chaine(4,virgule-5);
1981  Analyze(ctemp.Data(),err,offset); if (err) return;
1982  UInt_t leftopr = fNoper-1;
1983  ctemp = chaine(virgule,lchain-virgule-1);
1984  Analyze(ctemp.Data(),err,offset); if (err) return;
1985  fExpr[fNoper] = "min";
1986  actionCode = kmin;
1987  SetAction(fNoper,actionCode,actionParam);
1988  fNoper++;
1989  if (!CheckOperands(leftopr,fNoper-1,err)) return;
1990  }
1991  } else if (chaine(0,4) == "max(") {
1992  compt = 4; nomb = 0; virgule = 0; nest=0;
1993  while(compt != lchain) {
1994  compt++;
1995  if (chaine(compt-1,1) == "(") nest++;
1996  else if (chaine(compt-1,1) == ")") nest--;
1997  else if (chaine(compt-1,1) == "," && nest==0) {
1998  nomb++;
1999  if (nomb == 1 && virgule == 0) virgule = compt;
2000  }
2001  }
2002  if (nomb != 1) {
2003  err = 44; // There are plus or minus than 2 arguments for min
2004  err_hint = 3;
2005  }
2006  else {
2007  ctemp = chaine(4,virgule-5);
2008  Analyze(ctemp.Data(),err,offset); if (err) return;
2009  UInt_t leftopr = fNoper-1;
2010  ctemp = chaine(virgule,lchain-virgule-1);
2011  Analyze(ctemp.Data(),err,offset); if (err) return;
2012  fExpr[fNoper] = "max";
2013  actionCode = kmax;
2014  SetAction(fNoper,actionCode,actionParam);
2015  fNoper++;
2016  if (!CheckOperands(leftopr,fNoper-1,err)) return;
2017  }
2018 
2019  } else if (chaine(0,6) == "atan2(") {
2020  compt = 6; nomb = 0; virgule = 0; nest=0;
2021  while(compt != lchain) {
2022  compt++;
2023  if (chaine(compt-1,1) == "(") nest++;
2024  else if (chaine(compt-1,1) == ")") nest--;
2025  else if (chaine(compt-1,1) == "," && nest==0) {
2026  nomb++;
2027  if (nomb == 1 && virgule == 0) virgule = compt;
2028  }
2029  }
2030  if (nomb != 1) err = 21; //{ There are plus or minus than 2 arguments for atan2
2031  else {
2032  ctemp = chaine(6,virgule-7);
2033  Analyze(ctemp.Data(),err,offset); if (err) return;
2034  UInt_t leftopr = fNoper-1;
2035  ctemp = chaine(virgule,lchain-virgule-1);
2036  Analyze(ctemp.Data(),err,offset); if (err) return;
2037  fExpr[fNoper] = "atan2";
2038  actionCode = katan2;
2039  SetAction(fNoper,actionCode,actionParam);
2040  fNoper++;
2041  if (!CheckOperands(leftopr,fNoper-1,err)) return;
2042  }
2043  } else if (chaine(0,5) == "fmod(") {
2044  compt = 5; nomb = 0; virgule = 0; nest=0;
2045  while(compt != lchain) {
2046  compt++;
2047  if (chaine(compt-1,1) == "(") nest++;
2048  else if (chaine(compt-1,1) == ")") nest--;
2049  else if (chaine(compt-1,1) == "," && nest==0) {
2050  nomb++;
2051  if (nomb == 1 && virgule == 0) virgule = compt;
2052  }
2053  }
2054  if (nomb != 1) {
2055  err = 44; // There are plus or minus than 2 arguments for fmod
2056  err_hint = 4;
2057  }
2058  else {
2059  ctemp = chaine(5,virgule-6);
2060  Analyze(ctemp.Data(),err,offset); if (err) return;
2061  UInt_t leftopr = fNoper-1;
2062  ctemp = chaine(virgule,lchain-virgule-1);
2063  Analyze(ctemp.Data(),err,offset); if (err) return;
2064  fExpr[fNoper] = "fmod";
2065  actionCode = kfmod;
2066  SetAction(fNoper,actionCode,actionParam);
2067  fNoper++;
2068  if (!CheckOperands(leftopr,fNoper-1,err)) return;
2069  }
2070  } else if (AnalyzeFunction(chaine,err,offset) || err) { // The '||err' is to grab an error coming from AnalyzeFunction
2071  if (err) {
2072  chaine_error = chaine;
2073  } else {
2074  // We have a function call. Note that all the work was already,
2075  // eventually done in AnalyzeFunction
2076  //fprintf(stderr,"We found a foreign function in %s\n",chaine.Data());
2077  }
2078  } else if (chaine(0,1) == "[" && chaine(lchain-1,1) == "]") {
2079  fExpr[fNoper] = chaine;
2080  fNoper++;
2081  ctemp = chaine(1,lchain-2);
2082  for (j=0; j<ctemp.Length(); j++) {
2083  t=ctemp[j];
2084  if (strchr("0123456789",t)==0 && (ctemp(j,1)!="+" || j!=0)) {
2085  err=20;
2086  chaine_error=chaine1ST; // le numero ? de par[?] n'est pas un entier }
2087  }
2088  }
2089  if (!err) {
2090  sscanf(ctemp.Data(),"%d",&valeur);
2091  actionCode = kParameter;
2092  actionParam = offset + valeur;
2093  SetAction(fNoper-1, actionCode, actionParam);
2094  fExpr[fNoper-1] = "[";
2095  fExpr[fNoper-1] = (fExpr[fNoper-1] + (long int)(valeur+offset)) + "]";
2096  }
2097  } else if (chaine == "pi") {
2098  fExpr[fNoper] = "pi";
2099  actionCode = kpi;
2100  SetAction(fNoper,actionCode,actionParam);
2101  fNoper++;
2102  }
2103  else {
2104 
2105  // None of the above.
2106 
2107  err = 30;
2108  }
2109  }
2110  }
2111  }
2112  }
2113  }
2114  }
2115 
2116  // Overflows
2117  if (fNoper>=gMAXOP) err=6; // too many operators
2118 
2119  }
2120 
2121  // errors!
2122  if (err>1) {
2123  TString er = "";
2124  chaine_error = "\""+chaine_error+"\"";
2125  switch(err) {
2126  case 2 : er = " Invalid Floating Point Operation"; break;
2127  case 4 : er = " Empty String"; break;
2128  case 5 : er = " Invalid Syntax " + chaine_error; break;
2129  case 6 : er = " Too many operators !"; break;
2130  case 7 : er = " Too many parameters !"; break;
2131  case 10 : er = " z specified but not x and y"; break;
2132  case 11 : er = " z and y specified but not x"; break;
2133  case 12 : er = " y specified but not x"; break;
2134  case 13 : er = " z and x specified but not y"; break;
2135  case 20 : er = " Non integer value for parameter number : " + chaine_error; break;
2136  case 21 : er = " ATAN2 requires two arguments"; break;
2137  case 22 : er = " POW requires two arguments"; break;
2138  case 23 : er = " Degree of polynomial not specified"; break;
2139  case 24 : er = " Degree of polynomial must be positive"; break;
2140  case 25 : er = " Degree of polynomial must be less than 20"; break;
2141  case 26 : er = " Unknown name : " + chaine_error; break;
2142  case 27 : er = " Too many constants in expression"; break;
2143  case 28 : er = " strstr requires two arguments"; break;
2144  case 29 : er = " TFormula can only call interpreted and compiled functions that return a numerical type: " + chaine_error; break;
2145  case 30 : er = " Bad numerical expression : " + chaine_error; break;
2146  case 31 : er = " Part of the Variable " + chaine_error; er += " exists but some of it is not accessible or useable"; break;
2147  case 40 : er = " '(' is expected"; break;
2148  case 41 : er = " ')' is expected"; break;
2149  case 42 : er = " '[' is expected"; break;
2150  case 43 : er = " ']' is expected"; break;
2151  case 44 : er = " The function '" + chaine(0,err_hint) + "' requires two arguments."; break;
2152  case 45 : er = "The operator " + chaine_error + " requires a numerical operand."; break;
2153  case 46 : er = "Both operands of the operator " + chaine_error + " have to be either numbers or strings."; break;
2154  case 47 : er = chaine_error + " requires 2 string arguments"; break;
2155  }
2156  Error("Compile", "%s", er.Data());
2157  err=1;
2158  }
2159 
2160 }
2161 
2162 ////////////////////////////////////////////////////////////////////////////////
2163 /// Check whether the operand at 'oper-1' is compatible with the operation
2164 /// at 'oper'.
2165 
2166 Bool_t TFormula::CheckOperands(Int_t oper, Int_t &err)
2167 {
2168  if ( IsString(oper-1) && !StringToNumber(oper-1) ) {
2169  Error("Compile","\"%s\" requires a numerical operand.",fExpr[oper].Data());
2170  err = 45;
2171  return kFALSE;
2172  }
2173  return kTRUE;
2174 }
2175 
2176 ////////////////////////////////////////////////////////////////////////////////
2177 /// Check whether the operands at 'leftoper' and 'oper-1' are compatible with
2178 /// the operation at 'oper'.
2179 
2180 Bool_t TFormula::CheckOperands(Int_t leftoper, Int_t oper, Int_t &err)
2181 {
2182  if ( IsString(oper-1) || IsString(leftoper) ) {
2183  if (IsString(oper-1) && StringToNumber(oper-1)) {
2184  return kTRUE;
2185  }
2186  if (IsString(leftoper) && StringToNumber(leftoper)) {
2187  return kTRUE;
2188  }
2189  Error("Compile","\"%s\" requires two numerical operands.",fExpr[oper].Data());
2190  err = 46;
2191  return kFALSE;
2192  }
2193  return kTRUE;
2194 }
2195 
2196 ////////////////////////////////////////////////////////////////////////////////
2197 /// Try to 'demote' a string into an array bytes. If this is not possible,
2198 /// return false.
2199 
2200 Bool_t TFormula::StringToNumber(Int_t /* code */)
2201 {
2202  // In TFormula proper, we can not handle array of bytes ...
2203  return kFALSE;
2204 }
2205 
2206 ////////////////////////////////////////////////////////////////////////////////
2207 /// Resets the objects.
2208 ///
2209 /// Resets the object to its state before compilation.
2210 
2211 void TFormula::Clear(Option_t * /*option*/ )
2212 {
2213  ClearFormula();
2214 }
2215 
2216 ////////////////////////////////////////////////////////////////////////////////
2217 /// Resets the objects.
2218 ///
2219 /// Resets the object to its state before compilation.
2220 
2221 void TFormula::ClearFormula(Option_t * /*option*/ )
2222 {
2223  fNdim = 0;
2224  fNpar = 0;
2225  fNoper = 0;
2226  fNconst = 0;
2227  fNumber = 0;
2228  fNstring= 0;
2229  fNval = 0;
2230 
2231  if (fExpr) { delete [] fExpr; fExpr = 0;}
2232  if (fNames) { delete [] fNames; fNames = 0;}
2233  if (fOper) { delete [] fOper; fOper = 0;}
2234  if (fConst) { delete [] fConst; fConst = 0;}
2235  if (fParams) { delete [] fParams; fParams = 0;}
2236  fFunctions.Delete();
2237  fLinearParts.Delete();
2238  //
2239  //MI change
2240  if (fPredefined) { delete [] fPredefined; fPredefined = 0;}
2241  if (fOperOffset) { delete [] fOperOffset; fOperOffset = 0;}
2242  if (fExprOptimized) { delete [] fExprOptimized; fExprOptimized = 0;}
2243  if (fOperOptimized) { delete [] fOperOptimized; fOperOptimized = 0;}
2244  // should we also remove the object from the list?
2245  // gROOT->GetListOfFunctions()->Remove(this);
2246  // if we don't, what happens if it fails the new compilation?
2247 }
2248 
2249 namespace {
2250  template <class T>
2251  inline static void ResizeArrayIfAllocated(T*& oldArray, int newSize){
2252 
2253  // Don't do anything in this case.
2254  if (!oldArray || newSize <=0) return;
2255 
2256  T* newArray = new T[newSize];
2257  std::copy(oldArray, oldArray+newSize, newArray);
2258  delete [] oldArray;
2259  oldArray = newArray;
2260  }
2261 }
2262 
2263 ////////////////////////////////////////////////////////////////////////////////
2264 /// Compile expression already stored in fTitle.
2265 ///
2266 /// Loop on all subexpressions of formula stored in fTitle
2267 ///
2268 /// If you overload this member function, you also HAVE TO
2269 /// never call the constructor:
2270 ///
2271 /// ~~~ {.cpp}
2272 /// TFormula::TFormula(const char *name,const char *expression)
2273 /// ~~~
2274 ///
2275 /// and write your own constructor
2276 ///
2277 /// ~~~ {.cpp}
2278 /// MyClass::MyClass(const char *name,const char *expression) : TFormula()
2279 /// ~~~
2280 ///
2281 /// which has to call the TFormula default constructor and whose implementation
2282 /// should be similar to the implementation of the normal TFormula constructor
2283 ///
2284 /// This is necessary because the normal TFormula constructor call indirectly
2285 /// the virtual member functions Analyze, DefaultString, DefaultValue
2286 /// and DefaultVariable.
2287 ///
2288 /// \image html TFormula_compile.png
2289 
2290 Int_t TFormula::Compile(const char *expression)
2291 {
2292  Int_t i,j,lc,valeur,err;
2293  TString ctemp;
2294 
2295  ClearFormula();
2296 
2297  // If expression is not empty, take it, otherwise take the title
2298  if (strlen(expression)) SetTitle(expression);
2299 
2300  TString chaine = GetTitle();
2301 
2302  if (chaine.Contains(";")) {
2303  char *sctemp = new char[chaine.Length()+1];
2304  strlcpy(sctemp,chaine.Data(),chaine.Length()+1);
2305  char *semicol = (char*)strstr(sctemp,";");
2306  if (semicol) *semicol = 0;
2307  chaine = sctemp;
2308  delete [] sctemp;
2309  }
2310 
2311  // if the function is linear, process it and fill the array of linear parts
2312  if (TestBit(kLinear)){
2313  ProcessLinear(chaine);
2314  }
2315 
2316  // see static function SetMaxima to change the gMAX.. default values
2317  fExpr = new TString[gMAXOP];
2318  fConst = new Double_t[gMAXCONST];
2319  fParams = new Double_t[gMAXPAR];
2320  fNames = new TString[gMAXPAR];
2321  fOper = new Int_t[gMAXOP];
2322  for (i=0; i<gMAXPAR; i++) {
2323  fParams[i] = 0;
2324  fNames[i] = "";
2325  }
2326  for (i=0; i<gMAXOP; i++) {
2327  fExpr[i] = "";
2328  fOper[i] = 0;
2329  }
2330  for (i=0; i<gMAXCONST; i++)
2331  fConst[i] = 0;
2332 
2333  // Substitution of some operators to C++ style
2334  Bool_t inString = false;
2335  for (i=1; i<=chaine.Length(); i++) {
2336  lc =chaine.Length();
2337  if (chaine(i-1,1) == "\"") inString = !inString;
2338  if (inString) continue;
2339  if (chaine(i-1,2) == "**") {
2340  chaine = chaine(0,i-1) + "^" + chaine(i+1,lc-i-1);
2341  i=0;
2342  } else if (chaine(i-1,2) == "++") {
2343  chaine = chaine(0,i) + chaine(i+1,lc-i-1);
2344  i=0;
2345  } else if (chaine(i-1,2) == "+-" || chaine(i-1,2) == "-+") {
2346  chaine = chaine(0,i-1) + "-" + chaine(i+1,lc-i-1);
2347  i=0;
2348  } else if (chaine(i-1,2) == "--") {
2349  chaine = chaine(0,i-1) + "+" + chaine(i+1,lc-i-1);
2350  i=0;
2351  } else if (chaine(i-1,2) == "->") {
2352  chaine = chaine(0,i-1) + "." + chaine(i+1,lc-i-1);
2353  i=0;
2354  } else if (chaine(i-1,1) == "[") {
2355  for (j=1;j<=chaine.Length()-i;j++) {
2356  if (chaine(j+i-1,1) == "]" || j+i > chaine.Length()) break;
2357  }
2358  ctemp = chaine(i,j-1);
2359  valeur=0;
2360  sscanf(ctemp.Data(),"%d",&valeur);
2361  if (valeur >= fNpar) fNpar = valeur+1;
2362  } else if (chaine(i-1,1) == " ") {
2363  chaine = chaine(0,i-1)+chaine(i,lc-i);
2364  i=0;
2365  }
2366  }
2367  err = 0;
2368  Analyze((const char*)chaine,err);
2369 
2370  // if no parameters delete arrays fParams and fNames
2371  if (!fNpar) {
2372  delete [] fParams; fParams = 0;
2373  delete [] fNames; fNames = 0;
2374  }
2375 
2376  // if no errors, copy local parameters to formula objects
2377  if (!err) {
2378  if (fNdim <= 0) fNdim = 1;
2379  if (chaine.Length() > 4)
2380  {
2381  if ( GetNumber() != 400 &&
2382  GetNumber() != 410 &&
2383  GetNumber() != 110 )
2384  SetNumber(0);
2385  else if ( GetNumber() == 110 && chaine.Length() > 6 )
2386  SetNumber(0);
2387  else if ( GetNumber() == 410 && chaine.Length() > 8 )
2388  SetNumber(0);
2389  }
2390  // if formula is a gaussian, set parameter names
2391  if (GetNumber() == 100) {
2392  SetParName(0,"Constant");
2393  SetParName(1,"Mean");
2394  SetParName(2,"Sigma");
2395  }
2396  // if formula is a 2D gaussian, set parameter names
2397  if (GetNumber() == 110){
2398  SetParName(0,"Constant");
2399  SetParName(1,"MeanX");
2400  SetParName(2,"SigmaX");
2401  SetParName(3,"MeanY");
2402  SetParName(4,"SigmaY");
2403  }
2404  // if formula is an exponential, set parameter names
2405  if (GetNumber() == 200) {
2406  SetParName(0,"Constant");
2407  SetParName(1,"Slope");
2408  }
2409  // if formula is a polynom, set parameter names
2410  if (GetNumber() == 300+fNpar) {
2411  for (i=0;i<fNpar;i++) SetParName(i,Form("p%d",i));
2412  }
2413  // if formula is a landau, set parameter names
2414  if (GetNumber() == 400) {
2415  SetParName(0,"Constant");
2416  SetParName(1,"MPV");
2417  SetParName(2,"Sigma");
2418  }
2419  // if formula is a 2D landau, set parameter names
2420  if (GetNumber() == 410) {
2421  SetParName(0,"Constant");
2422  SetParName(1,"MPVX");
2423  SetParName(2,"SigmaX");
2424  SetParName(3,"MPVY");
2425  SetParName(4,"SigmaY");
2426  }
2427  }
2428 
2429  // Here we shrink the arrays allocated like this:
2430  // fExpr = new TString[gMAXOP];
2431  // fConst = new Double_t[gMAXCONST];
2432  // fParams = new Double_t[gMAXPAR];
2433  // fNames = new TString[gMAXPAR];
2434  // fOper = new Int_t[gMAXOP];
2435  // fParams and fNames may be already 0, so we have to check.
2436  if (!err){
2437  ResizeArrayIfAllocated(fExpr, fNoper);
2438  ResizeArrayIfAllocated(fConst, fNconst);
2439  ResizeArrayIfAllocated(fParams, fNpar);
2440  ResizeArrayIfAllocated(fNames, fNpar);
2441  ResizeArrayIfAllocated(fOper, fNoper);
2442  }
2443 
2444 
2445  if (err) { fNdim = 0; return 1; }
2446  // Convert(5);
2447  //
2448  //MI change
2449  if (!IsA()->GetBaseClass("TTreeFormula")) {
2450  Optimize();
2451  }
2452  //
2453  return 0;
2454 }
2455 
2456 ////////////////////////////////////////////////////////////////////////////////
2457 /// Copy this formula.
2458 
2459 void TFormula::Copy(TObject &obj) const
2460 {
2461  Int_t i;
2462  ((TFormula&)obj).ClearFormula();
2463  TNamed::Copy(obj);
2464  ((TFormula&)obj).fNdim = fNdim;
2465  ((TFormula&)obj).fNpar = fNpar;
2466  ((TFormula&)obj).fNoper = fNoper;
2467  ((TFormula&)obj).fNconst = fNconst;
2468  ((TFormula&)obj).fNumber = fNumber;
2469  ((TFormula&)obj).fNval = fNval;
2470  ((TFormula&)obj).fExpr = 0;
2471  ((TFormula&)obj).fConst = 0;
2472  ((TFormula&)obj).fParams = 0;
2473  ((TFormula&)obj).fNames = 0;
2474  if (fExpr && fNoper) {
2475  ((TFormula&)obj).fExpr = new TString[fNoper];
2476  for (i=0;i<fNoper;i++) ((TFormula&)obj).fExpr[i] = fExpr[i];
2477  }
2478  if (fOper && fNoper) {
2479  ((TFormula&)obj).fOper = new Int_t[fNoper];
2480  for (i=0;i<fNoper;i++) ((TFormula&)obj).fOper[i] = fOper[i];
2481  }
2482  if (fConst && fNconst) {
2483  ((TFormula&)obj).fConst = new Double_t[fNconst];
2484  for (i=0;i<fNconst;i++) ((TFormula&)obj).fConst[i] = fConst[i];
2485  }
2486  if (fParams && fNpar) {
2487  ((TFormula&)obj).fParams = new Double_t[fNpar];
2488  for (i=0;i<fNpar;i++) ((TFormula&)obj).fParams[i] = fParams[i];
2489  }
2490  if (fNames && fNpar) {
2491  ((TFormula&)obj).fNames = new TString[fNpar];
2492  for (i=0;i<fNpar;i++) ((TFormula&)obj).fNames[i] = fNames[i];
2493  }
2494 
2495  TIter next(&fFunctions);
2496  TObject *fobj;
2497  while ( (fobj = next()) ) {
2498  ((TFormula&)obj).fFunctions.Add( fobj->Clone() );
2499  }
2500  //
2501  // MI change
2502  //
2503  //
2504  if (fNoper) {
2505  if(fExprOptimized) {
2506  ((TFormula&)obj).fExprOptimized = new TString[fNoper];
2507  for (i=0;i<fNoper;i++) ((TFormula&)obj).fExprOptimized[i] = fExprOptimized[i];
2508  }
2509  if (fOperOptimized) {
2510  ((TFormula&)obj).fOperOptimized = new Int_t[fNoper];
2511  for (i=0;i<fNoper;i++) ((TFormula&)obj).fOperOptimized[i] = fOperOptimized[i];
2512  }
2513  if (fPredefined) {
2514  ((TFormula&)obj).fPredefined = new ROOT::v5::TFormulaPrimitive*[fNoper];
2515  for (i=0;i<fNoper;i++) {((TFormula&)obj).fPredefined[i] = fPredefined[i];}
2516  }
2517  if (fOperOffset) {
2518  ((TFormula&)obj).fOperOffset = new TOperOffset[fNoper];
2519  for (i=0;i<fNoper;i++) {((TFormula&)obj).fOperOffset[i] = fOperOffset[i];}
2520  }
2521  }
2522  ((TFormula&)obj).fNOperOptimized = fNOperOptimized;
2523  ((TFormula&)obj).fOptimal = fOptimal;
2524 
2525 }
2526 
2527 ////////////////////////////////////////////////////////////////////////////////
2528 /// Return address of string corresponding to special code.
2529 ///
2530 /// This member function is inactive in the TFormula class.
2531 /// It may be redefined in derived classes.
2532 ///
2533 /// If you overload this member function, you also HAVE TO
2534 /// never call the constructor:
2535 ///
2536 /// ~~~ {.cpp}
2537 /// TFormula::TFormula(const char *name,const char *expression)
2538 /// ~~~
2539 ///
2540 /// and write your own constructor
2541 ///
2542 /// ~~~ {.cpp}
2543 /// MyClass::MyClass(const char *name,const char *expression) : TFormula()
2544 /// ~~~
2545 ///
2546 /// which has to call the TFormula default constructor and whose implementation
2547 /// should be similar to the implementation of the normal TFormula constructor
2548 ///
2549 /// This is necessary because the normal TFormula constructor call indirectly
2550 /// the virtual member functions Analyze, DefaultString, DefaultValue
2551 /// and DefaultVariable.
2552 
2553 char *TFormula::DefinedString(Int_t)
2554 {
2555  return 0;
2556 }
2557 
2558 ////////////////////////////////////////////////////////////////////////////////
2559 /// Return value corresponding to special code.
2560 ///
2561 /// This member function is inactive in the TFormula class.
2562 /// It may be redefined in derived classes.
2563 ///
2564 /// If you overload this member function, you also HAVE TO
2565 /// never call the constructor:
2566 ///
2567 /// ~~~ {.cpp}
2568 /// TFormula::TFormula(const char *name,const char *expression)
2569 /// ~~~
2570 ///
2571 /// and write your own constructor
2572 ///
2573 /// ~~~ {.cpp}
2574 /// MyClass::MyClass(const char *name,const char *expression) : TFormula()
2575 /// ~~~
2576 ///
2577 /// which has to call the TFormula default constructor and whose implementation
2578 /// should be similar to the implementation of the normal TFormula constructor
2579 ///
2580 /// This is necessary because the normal TFormula constructor call indirectly
2581 /// the virtual member functions Analyze, DefaultString, DefaultValue
2582 /// and DefaultVariable.
2583 
2584 Double_t TFormula::DefinedValue(Int_t)
2585 {
2586  return 0;
2587 }
2588 
2589 ////////////////////////////////////////////////////////////////////////////////
2590 /// Check if expression is in the list of defined variables.
2591 ///
2592 /// This member function can be overloaded in derived classes
2593 ///
2594 /// If you overload this member function, you also HAVE TO
2595 /// never call the constructor:
2596 ///
2597 /// ~~~ {.cpp}
2598 /// TFormula::TFormula(const char *name,const char *expression)
2599 /// ~~~
2600 ///
2601 /// and write your own constructor
2602 ///
2603 /// ~~~ {.cpp}
2604 /// MyClass::MyClass(const char *name,const char *expression) : TFormula()
2605 /// ~~~
2606 ///
2607 /// which has to call the TFormula default constructor and whose implementation
2608 /// should be similar to the implementation of the normal TFormula constructor
2609 ///
2610 /// This is necessary because the normal TFormula constructor call indirectly
2611 /// the virtual member functions Analyze, DefaultString, DefaultValue
2612 /// and DefaultVariable.
2613 ///
2614 /// The expected returns values are
2615 /// - -2 : the name has been recognized but won't be usable
2616 /// - -1 : the name has not been recognized
2617 /// - >=0 : the name has been recognized, return the action parameter.
2618 
2619 Int_t TFormula::DefinedVariable(TString &chaine,Int_t &action)
2620 {
2621  action = kVariable;
2622  if (chaine == "x") {
2623  if (fNdim < 1) fNdim = 1;
2624  return 0;
2625  } else if (chaine == "y") {
2626  if (fNdim < 2) fNdim = 2;
2627  return 1;
2628  } else if (chaine == "z") {
2629  if (fNdim < 3) fNdim = 3;
2630  return 2;
2631  } else if (chaine == "t") {
2632  if (fNdim < 4) fNdim = 4;
2633  return 3;
2634  }
2635  // MI change
2636  // extended defined variable (MI)
2637  //
2638  if (chaine.Data()[0]=='x'){
2639  if (chaine.Data()[1]=='[' && chaine.Data()[3]==']'){
2640  const char ch0 = '0';
2641  Int_t dim = chaine.Data()[2]-ch0;
2642  if (dim<0) return -1;
2643  if (dim>9) return -1;
2644  if (fNdim<=dim) fNdim = dim+1;
2645  return dim;
2646  }
2647  if (chaine.Data()[1]=='[' && chaine.Data()[4]==']'){
2648  const char ch0 = '0';
2649  Int_t dim = (chaine.Data()[2]-ch0)*10+(chaine.Data()[3]-ch0);
2650  if (dim<0) return -1;
2651  if (dim>99) return -1;
2652  if (fNdim<=dim) fNdim = dim+1;
2653  return dim;
2654  }
2655  }
2656  return -1;
2657 }
2658 
2659 ////////////////////////////////////////////////////////////////////////////////
2660 /// Evaluate this formula.
2661 ///
2662 /// The current value of variables x,y,z,t is passed through x, y, z and t.
2663 /// The parameters used will be the ones in the array params if params is given
2664 /// otherwise parameters will be taken from the stored data members fParams
2665 
2666 Double_t TFormula::Eval(Double_t x, Double_t y, Double_t z, Double_t t) const
2667 {
2668  Double_t xx[4];
2669  xx[0] = x;
2670  xx[1] = y;
2671  xx[2] = z;
2672  xx[3] = t;
2673  return ((TFormula*)this)->EvalPar(xx);
2674 }
2675 
2676 ////////////////////////////////////////////////////////////////////////////////
2677 /// Evaluate this formula.
2678 ///
2679 /// The current value of variables x,y,z,t is passed through the pointer x.
2680 /// The parameters used will be the ones in the array params if params is given
2681 /// otherwise parameters will be taken from the stored data members fParams
2682 ///
2683 /// \image html TFormula_eval.png
2684 
2685 Double_t TFormula::EvalParOld(const Double_t *x, const Double_t *uparams)
2686 {
2687  Int_t i,j;
2688  // coverity[uninit] the tab value of tab is guaranteed to be set properly by the control flow.
2689  Double_t tab[kMAXFOUND];
2690  const char *stringStack[gMAXSTRINGFOUND] = {0};
2691  Double_t param_calc[kMAXFOUND];
2692  char *string_calc[gMAXSTRINGFOUND] = {0};
2693  Int_t precalculated = 0;
2694  Int_t precalculated_str = 0;
2695  Double_t *params;
2696 
2697  if (uparams) {
2698  params = const_cast<Double_t*>(uparams);
2699  } else {
2700  params = fParams;
2701  }
2702  UInt_t pos = 0;
2703  UInt_t strpos = 0;
2704 
2705  for (i=0; i<fNoper; ++i) {
2706 
2707  const int oper = fOper[i];
2708  const int opcode = oper >> kTFOperShift;
2709 
2710  switch(opcode) {
2711 
2712  case kParameter : { pos++; tab[pos-1] = params[ oper & kTFOperMask ]; continue; }
2713  case kConstant : { pos++; tab[pos-1] = fConst[ oper & kTFOperMask ]; continue; }
2714  case kVariable : { pos++; tab[pos-1] = x[ oper & kTFOperMask ]; continue; }
2715  case kStringConst: { strpos++; stringStack[strpos-1] = (char*)fExpr[i].Data(); pos++; tab[pos-1] = 0; continue; }
2716 
2717  case kAdd : pos--; tab[pos-1] += tab[pos]; continue;
2718  case kSubstract : pos--; tab[pos-1] -= tab[pos]; continue;
2719  case kMultiply : pos--; tab[pos-1] *= tab[pos]; continue;
2720  case kDivide : pos--; if (tab[pos] == 0) tab[pos-1] = 0; // division by 0
2721  else tab[pos-1] /= tab[pos];
2722  continue;
2723  case kModulo : {pos--;
2724  Long64_t int1((Long64_t)tab[pos-1]);
2725  Long64_t int2((Long64_t)tab[pos]);
2726  tab[pos-1] = Double_t(int1%int2);
2727  continue;}
2728 
2729  case kcos : tab[pos-1] = TMath::Cos(tab[pos-1]); continue;
2730  case ksin : tab[pos-1] = TMath::Sin(tab[pos-1]); continue;
2731  case ktan : if (TMath::Cos(tab[pos-1]) == 0) {tab[pos-1] = 0;} // { tangente indeterminee }
2732  else tab[pos-1] = TMath::Tan(tab[pos-1]);
2733  continue;
2734  case kacos : if (TMath::Abs(tab[pos-1]) > 1) {tab[pos-1] = 0;} // indetermination
2735  else tab[pos-1] = TMath::ACos(tab[pos-1]);
2736  continue;
2737  case kasin : if (TMath::Abs(tab[pos-1]) > 1) {tab[pos-1] = 0;} // indetermination
2738  else tab[pos-1] = TMath::ASin(tab[pos-1]);
2739  continue;
2740  case katan : tab[pos-1] = TMath::ATan(tab[pos-1]); continue;
2741  case kcosh : tab[pos-1] = TMath::CosH(tab[pos-1]); continue;
2742  case ksinh : tab[pos-1] = TMath::SinH(tab[pos-1]); continue;
2743  case ktanh : if (TMath::CosH(tab[pos-1]) == 0) {tab[pos-1] = 0;} // { tangente indeterminee }
2744  else tab[pos-1] = TMath::TanH(tab[pos-1]);
2745  continue;
2746  case kacosh: if (tab[pos-1] < 1) {tab[pos-1] = 0;} // indetermination
2747  else tab[pos-1] = TMath::ACosH(tab[pos-1]);
2748  continue;
2749  case kasinh: tab[pos-1] = TMath::ASinH(tab[pos-1]); continue;
2750  case katanh: if (TMath::Abs(tab[pos-1]) > 1) {tab[pos-1] = 0;} // indetermination
2751  else tab[pos-1] = TMath::ATanH(tab[pos-1]);
2752  continue;
2753  case katan2: pos--; tab[pos-1] = TMath::ATan2(tab[pos-1],tab[pos]); continue;
2754 
2755  case kfmod : pos--; tab[pos-1] = fmod(tab[pos-1],tab[pos]); continue;
2756  case kpow : pos--; tab[pos-1] = TMath::Power(tab[pos-1],tab[pos]); continue;
2757  case ksq : tab[pos-1] = tab[pos-1]*tab[pos-1]; continue;
2758  case ksqrt : tab[pos-1] = TMath::Sqrt(TMath::Abs(tab[pos-1])); continue;
2759 
2760  case kstrstr : strpos -= 2; pos-=2; pos++;
2761  if (strstr(stringStack[strpos],stringStack[strpos+1])) tab[pos-1]=1;
2762  else tab[pos-1]=0;
2763  continue;
2764 
2765  case kmin : pos--; tab[pos-1] = TMath::Min(tab[pos-1],tab[pos]); continue;
2766  case kmax : pos--; tab[pos-1] = TMath::Max(tab[pos-1],tab[pos]); continue;
2767 
2768  case klog : if (tab[pos-1] > 0) tab[pos-1] = TMath::Log(tab[pos-1]);
2769  else {tab[pos-1] = 0;} //{indetermination }
2770  continue;
2771  case kexp : { Double_t dexp = tab[pos-1];
2772  if (dexp < -700) {tab[pos-1] = 0; continue;}
2773  if (dexp > 700) {tab[pos-1] = TMath::Exp(700); continue;}
2774  tab[pos-1] = TMath::Exp(dexp); continue; }
2775  case klog10: if (tab[pos-1] > 0) tab[pos-1] = TMath::Log10(tab[pos-1]);
2776  else {tab[pos-1] = 0;} //{indetermination }
2777  continue;
2778 
2779  case kpi : pos++; tab[pos-1] = TMath::ACos(-1); continue;
2780 
2781  case kabs : tab[pos-1] = TMath::Abs(tab[pos-1]); continue;
2782  case ksign : if (tab[pos-1] < 0) tab[pos-1] = -1; else tab[pos-1] = 1; continue;
2783  case kint : tab[pos-1] = Double_t(Int_t(tab[pos-1])); continue;
2784 
2785  case kSignInv: tab[pos-1] = -1 * tab[pos-1]; continue;
2786 
2787  case krndm : pos++; tab[pos-1] = gRandom->Rndm(); continue;
2788 
2789  case kAnd : pos--; if (tab[pos-1]!=0 && tab[pos]!=0) tab[pos-1]=1;
2790  else tab[pos-1]=0;
2791  continue;
2792  case kOr : pos--; if (tab[pos-1]!=0 || tab[pos]!=0) tab[pos-1]=1;
2793  else tab[pos-1]=0;
2794  continue;
2795  case kEqual: pos--; if (tab[pos-1] == tab[pos]) tab[pos-1]=1;
2796  else tab[pos-1]=0;
2797  continue;
2798  case kNotEqual : pos--; if (tab[pos-1] != tab[pos]) tab[pos-1]=1;
2799  else tab[pos-1]=0;
2800  continue;
2801  case kLess : pos--; if (tab[pos-1] < tab[pos]) tab[pos-1]=1;
2802  else tab[pos-1]=0;
2803  continue;
2804  case kGreater : pos--; if (tab[pos-1] > tab[pos]) tab[pos-1]=1;
2805  else tab[pos-1]=0;
2806  continue;
2807 
2808  case kLessThan: pos--; if (tab[pos-1]<=tab[pos]) tab[pos-1]=1;
2809  else tab[pos-1]=0;
2810  continue;
2811  case kGreaterThan: pos--; if (tab[pos-1]>=tab[pos]) tab[pos-1]=1;
2812  else tab[pos-1]=0;
2813  continue;
2814  case kNot : if (tab[pos-1]!=0) tab[pos-1] = 0; else tab[pos-1] = 1;
2815  continue;
2816 
2817  case kStringEqual : strpos -= 2; pos -=2 ; pos++;
2818  if (!strcmp(stringStack[strpos+1],stringStack[strpos])) tab[pos-1]=1;
2819  else tab[pos-1]=0;
2820  continue;
2821  case kStringNotEqual: strpos -= 2; pos -= 2; pos++;
2822  if (strcmp(stringStack[strpos+1],stringStack[strpos])) tab[pos-1]=1;
2823  else tab[pos-1]=0;
2824  continue;
2825 
2826  case kBitAnd : pos--; tab[pos-1]= ((Int_t) tab[pos-1]) & ((Int_t) tab[pos]); continue;
2827  case kBitOr : pos--; tab[pos-1]= ((Int_t) tab[pos-1]) | ((Int_t) tab[pos]); continue;
2828  case kLeftShift : pos--; tab[pos-1]= ((Int_t) tab[pos-1]) <<((Int_t) tab[pos]); continue;
2829  case kRightShift: pos--; tab[pos-1]= ((Int_t) tab[pos-1]) >>((Int_t) tab[pos]); continue;
2830 
2831  case kJump : i = (oper & kTFOperMask); continue;
2832  case kJumpIf : pos--; if (!tab[pos]) i = (oper & kTFOperMask); continue;
2833 
2834  case kBoolOptimize: {
2835  // boolean operation optimizer
2836 
2837  int param = (oper & kTFOperMask);
2838  Bool_t skip = kFALSE;
2839  int op = param % 10; // 1 is && , 2 is ||
2840 
2841  if (op == 1 && (!tab[pos-1]) ) {
2842  // &&: skip the right part if the left part is already false
2843 
2844  skip = kTRUE;
2845 
2846  // Preserve the existing behavior (i.e. the result of a&&b is
2847  // either 0 or 1)
2848  tab[pos-1] = 0;
2849 
2850  } else if (op == 2 && tab[pos-1] ) {
2851  // ||: skip the right part if the left part is already true
2852 
2853  skip = kTRUE;
2854 
2855  // Preserve the existing behavior (i.e. the result of a||b is
2856  // either 0 or 1)
2857  tab[pos-1] = 1;
2858 
2859  }
2860 
2861  if (skip) {
2862  int toskip = param / 10;
2863  i += toskip;
2864  }
2865  continue;
2866  }
2867 
2868  }
2869 
2870  switch(opcode) {
2871 
2872  #define R__EXPO(var) \
2873  { \
2874  pos++; int param = (oper & kTFOperMask); \
2875  tab[pos-1] = TMath::Exp(params[param]+params[param+1]*x[var]); \
2876  continue; \
2877  }
2878  // case kexpo:
2879  case kxexpo: R__EXPO(0);
2880  case kyexpo: R__EXPO(1);
2881  case kzexpo: R__EXPO(2);
2882  case kxyexpo:{ pos++; int param = (oper & kTFOperMask);
2883  tab[pos-1] = TMath::Exp(params[param]+params[param+1]*x[0]+params[param+2]*x[1]);
2884  continue; }
2885 
2886  #define R__GAUS(var) \
2887  { \
2888  pos++; int param = (oper & kTFOperMask); \
2889  tab[pos-1] = params[param]*TMath::Gaus(x[var],params[param+1],params[param+2],IsNormalized()); \
2890  continue; \
2891  }
2892 
2893  // case kgaus:
2894  case kxgaus: R__GAUS(0);
2895  case kygaus: R__GAUS(1);
2896  case kzgaus: R__GAUS(2);
2897  case kxygaus: {pos++; int param = (oper & kTFOperMask);
2898  Double_t intermede1;
2899  if (params[param+2] == 0) {
2900  intermede1=1e10;
2901  } else {
2902  intermede1=Double_t((x[0]-params[param+1])/params[param+2]);
2903  }
2904  Double_t intermede2;
2905  if (params[param+4] == 0) {
2906  intermede2=1e10;
2907  } else {
2908  intermede2=Double_t((x[1]-params[param+3])/params[param+4]);
2909  }
2910  tab[pos-1] = params[param]*TMath::Exp(-0.5*(intermede1*intermede1+intermede2*intermede2));
2911  continue; }
2912 
2913  #define R__LANDAU(var) \
2914  { \
2915  pos++; const int param = (oper & kTFOperMask); \
2916  tab[pos-1] = params[param]*TMath::Landau(x[var],params[param+1],params[param+2],IsNormalized()); \
2917  continue; \
2918  }
2919  // case klandau:
2920  case kxlandau: R__LANDAU(0);
2921  case kylandau: R__LANDAU(1);
2922  case kzlandau: R__LANDAU(2);
2923  case kxylandau: { pos++; int param = oper&0x7fffff /* ActionParams[i] */ ;
2924  Double_t intermede1=TMath::Landau(x[0], params[param+1], params[param+2],IsNormalized());
2925  Double_t intermede2=TMath::Landau(x[1], params[param+3], params[param+4],IsNormalized());
2926  tab[pos-1] = params[param]*intermede1*intermede2;
2927  continue;
2928  }
2929 
2930  #define R__POLY(var) \
2931  { \
2932  pos++; int param = (oper & kTFOperMask); \
2933  tab[pos-1] = 0; Double_t intermede = 1; \
2934  Int_t inter = param/100; /* arrondit */ \
2935  Int_t int1= param-inter*100-1; /* aucune simplification ! (sic) */ \
2936  for (j=0 ;j<inter+1;j++) { \
2937  tab[pos-1] += intermede*params[j+int1]; \
2938  intermede *= x[var]; \
2939  } \
2940  continue; \
2941  }
2942  // case kpol:
2943  case kxpol: R__POLY(0);
2944  case kypol: R__POLY(1);
2945  case kzpol: R__POLY(2);
2946 
2947  case kDefinedVariable : {
2948  if (!precalculated) {
2949  precalculated = 1;
2950  for(j=0;j<fNval;j++) param_calc[j]=DefinedValue(j);
2951  }
2952  pos++; tab[pos-1] = param_calc[(oper & kTFOperMask)];
2953  continue;
2954  }
2955 
2956  case kDefinedString : {
2957  int param = (oper & kTFOperMask);
2958  if (!precalculated_str) {
2959  precalculated_str=1;
2960  for (j=0;j<fNstring;j++) string_calc[j]=DefinedString(j);
2961  }
2962  strpos++; stringStack[strpos-1] = string_calc[param];
2963  pos++; tab[pos-1] = 0;
2964  continue;
2965  }
2966 
2967  case kFunctionCall: {
2968  // an external function call
2969 
2970  int param = (oper & kTFOperMask);
2971  int fno = param / 1000;
2972  int nargs = param % 1000;
2973 
2974  // Retrieve the function
2975  TMethodCall *method = (TMethodCall*)fFunctions.At(fno);
2976 
2977  // Set the arguments
2978  method->ResetParam();
2979  if (nargs) {
2980  UInt_t argloc = pos-nargs;
2981  for(j=0;j<nargs;j++,argloc++,pos--) {
2982  method->SetParam(tab[argloc]);
2983  }
2984  }
2985  pos++;
2986  Double_t ret;
2987  method->Execute(ret);
2988  tab[pos-1] = ret; // check for the correct conversion!
2989 
2990  continue;
2991  };
2992  }
2993  if (!TestBit(kOptimizationError)) {
2994  SetBit(kOptimizationError);
2995  Warning("EvalParOld","Found an unsupported opcode (%d)",oper >> kTFOperShift);
2996  }
2997  }
2998  Double_t result0 = tab[0];
2999  return result0;
3000 
3001 }
3002 
3003 ////////////////////////////////////////////////////////////////////////////////
3004 /// Reconstruct the formula expression from the internal TFormula member variables
3005 ///
3006 /// This function uses the internal member variables of TFormula to
3007 /// construct the mathematical expression associated with the TFormula
3008 /// instance. This function can be used to get an expanded version of the
3009 /// expression originally assigned to the TFormula instance, i.e. that
3010 /// the string returned by GetExpFormula() doesn't depend on other
3011 /// TFormula object names.
3012 ///
3013 /// if option contains "p" the returned string will contain the formula
3014 /// expression with symbolic parameters, eg [0] replaced by the actual value
3015 /// of the parameter. Example:
3016 /// if expression in formula is: "[0]*(x>-[1])+[2]*exp(-[3]*x)"
3017 /// and parameters are 3.25,-4.01,4.44,-0.04, GetExpFormula("p") will return:
3018 /// "(3.25*(x>+4.01))+(4.44*exp(+0.04*x))"
3019 
3020 TString TFormula::GetExpFormula(Option_t *option) const
3021 {
3022  if (fNoper>0) {
3023  TString* tab=new TString[fNoper];
3024  Bool_t* ismulti=new Bool_t[fNoper];
3025  Int_t spos=0;
3026 
3027  ismulti[0]=kFALSE;
3028  Int_t optype;
3029  Int_t j;
3030  Int_t ternaryend = -1;
3031  for(Int_t i=0;i<fNoper;i++){
3032  optype = GetAction(i);
3033 
3034  if (ternaryend==i) {
3035  // The ? and : have been added to tab[spos-2]
3036  if(ismulti[spos-1]){
3037  tab[spos-2]=tab[spos-2]+"("+tab[spos-1]+")";
3038  } else {
3039  tab[spos-2]=tab[spos-2]+tab[spos-1];
3040  }
3041  spos--;
3042  // Do not call continue since we need to
3043  // do the rest of the loop.
3044  }
3045 
3046  // Boolean optimization breakpoint
3047  if (optype==kBoolOptimize) { // -3) {
3048  continue;
3049  }
3050 
3051  //Sign inversion
3052  if (optype==kSignInv) { // -1) {
3053  tab[spos-1]="-("+tab[spos-1]+")";
3054  // i++;
3055  continue;
3056  }
3057 
3058  //Simple name (parameter,pol0,landau, etc)
3059  if (kexpo<=optype && optype<=kzpol) { // >=0) {
3060  tab[spos]=fExpr[i];
3061  ismulti[spos]=kFALSE;
3062  spos++;
3063  continue;
3064  }
3065 
3066  //constants, variables x,y,z,t, pi
3067  if ((optype<=151 && optype>=140 && optype!=145) || (optype == 40)) {
3068  tab[spos]=fExpr[i];
3069  ismulti[spos]=kFALSE;
3070  spos++;
3071  continue;
3072  }
3073 
3074  //Basic operators (+,-,*,/,==,^,etc)
3075  if(((optype>0 && optype<6) || optype==20 ||
3076  (((optype>59 && optype<69) || (optype >75 && optype<82)) && spos>=2))) {
3077  // if(optype==-20 && spos>=2){
3078  if(ismulti[spos-2]){
3079  tab[spos-2]="("+tab[spos-2]+")";
3080  }
3081  if(ismulti[spos-1]){
3082  tab[spos-2]+=fExpr[i]+("("+tab[spos-1]+")");
3083  }else{
3084  tab[spos-2]+=fExpr[i]+tab[spos-1];
3085  }
3086  ismulti[spos-2]=kTRUE;
3087  spos--;
3088  continue;
3089  }
3090  //Ternary condition
3091  if (optype==kJumpIf) {
3092  if(ismulti[spos-1]){
3093  tab[spos-1]="("+tab[spos-1]+")?";
3094  } else {
3095  tab[spos-1]=tab[spos-1]+"?";
3096  }
3097  continue;
3098  }
3099  if (optype==kJump) {
3100  if(ismulti[spos-1]){
3101  tab[spos-2]=tab[spos-2]+"("+tab[spos-1]+"):";
3102  } else {
3103  tab[spos-2]=tab[spos-2]+tab[spos-1]+":";
3104  }
3105  ternaryend = GetActionParam(i) + 1;
3106  spos--;
3107  continue;
3108  }
3109 
3110  //Functions
3111  int offset = 0;
3112  TString funcname = fExpr[i];
3113  if((optype>9 && optype<16) ||
3114  (optype>20 && optype<23) ||
3115  (optype>29 && optype<34) ||
3116  (optype>40 && optype<44) ||
3117  (optype>69 && optype<76) ||
3118  (optype==145)) {
3119  //Functions with the format func(x)
3120  offset = -1;
3121  }
3122 
3123  if((optype>15 && optype<20) ||
3124  (optype>22 && optype<26)) {
3125  //Functions with the format func(x,y)
3126  offset = -2;
3127  }
3128  if(optype==145) {
3129  int param = (fOper[i] & kTFOperMask);
3130  //int fno = param / 1000;
3131  int nargs = param % 1000;
3132  offset = -nargs;
3133  // The function name contains return type and parameters types we need
3134  // to trim them.
3135  int depth;
3136  for(j=0, depth=0;j<funcname.Length();++j) {
3137  switch (funcname[j]) {
3138  case '<':
3139  ++depth; break;
3140  case '>':
3141  --depth; break;
3142  case ' ':
3143  if (depth==0) {
3144  funcname.Remove(0,j+1);
3145  j = funcname.Length();
3146  break;
3147  }
3148  }
3149  }
3150  Ssiz_t ind = funcname.First('(');
3151  funcname.Remove(ind);
3152  }
3153  if (offset > 0) {
3154  Error("GetExpFormula","Internal error, number of argument found is %d",-offset);
3155  } else if (offset == 0) {
3156  tab[spos]=funcname+"()";
3157  ismulti[spos]=kFALSE;
3158  spos += 1;
3159  continue;
3160  } else if (offset<=0 && (spos+offset>=0)) {
3161  tab[spos+offset]=funcname+("("+tab[spos+offset]);
3162  for (j=offset+1; j<0; j++){
3163  tab[spos+offset]+=","+tab[spos+j];
3164  }
3165  tab[spos+offset]+=")";
3166  ismulti[spos+offset]=kFALSE;
3167  spos += offset+1;
3168  continue;
3169  }
3170  }
3171  if (ternaryend==fNoper) {
3172  // The ? and : have been added to tab[spos-2]
3173  if(ismulti[spos-1]){
3174  tab[spos-2]=tab[spos-2]+"("+tab[spos-1]+")";
3175  } else {
3176  tab[spos-2]=tab[spos-2]+tab[spos-1];
3177  }
3178  spos--;
3179  }
3180 
3181  TString ret = "";
3182  if (spos > 0) ret = tab[spos-1];
3183  delete[] tab;
3184  delete[] ismulti;
3185 
3186  //if option "p" is specified, return the real values of parameters instead of [0]
3187  TString opt = option;
3188  opt.ToLower();
3189  if (opt.Contains("p")) {
3190  char pb[13];
3191  char pbv[100];
3192  for (j=0;j<fNpar;j++) {
3193  snprintf(pb,sizeof(pb),"[%d]",j);
3194  snprintf(pbv,100,"%g",fParams[j]);
3195  ret.ReplaceAll(pb,pbv);
3196  }
3197  ret.ReplaceAll("--","+");
3198  ret.ReplaceAll("+-","-");
3199  }
3200  return ret;
3201  } else{
3202  TString ret="";
3203  return ret;
3204  }
3205 }
3206 
3207 ////////////////////////////////////////////////////////////////////////////////
3208 /// Return linear part.
3209 
3210 const TObject* TFormula::GetLinearPart(Int_t i)
3211 {
3212  if (!fLinearParts.IsEmpty())
3213  return fLinearParts.UncheckedAt(i);
3214  return 0;
3215 }
3216 
3217 ////////////////////////////////////////////////////////////////////////////////
3218 /// Return value of parameter number ipar.
3219 
3220 Double_t TFormula::GetParameter(Int_t ipar) const
3221 {
3222  if (ipar <0 || ipar >= fNpar) return 0;
3223  return fParams[ipar];
3224 }
3225 
3226 ////////////////////////////////////////////////////////////////////////////////
3227 /// Return value of parameter named parName.
3228 
3229 Double_t TFormula::GetParameter(const char *parName) const
3230 {
3231  const Double_t kNaN = 1e-300;
3232  Int_t index = GetParNumber(parName);
3233  if (index==-1) {
3234  Error("TFormula", "Parameter %s not found", parName);
3235  return kNaN;
3236  }
3237  return GetParameter(index);
3238 }
3239 
3240 ////////////////////////////////////////////////////////////////////////////////
3241 /// Return name of one parameter.
3242 
3243 const char *TFormula::GetParName(Int_t ipar) const
3244 {
3245  if (ipar <0 || ipar >= fNpar) return "";
3246  if (fNames[ipar].Length() > 0) return (const char*)fNames[ipar];
3247  return Form("p%d",ipar);
3248 }
3249 
3250 ////////////////////////////////////////////////////////////////////////////////
3251 /// Return parameter number by name.
3252 
3253 Int_t TFormula::GetParNumber(const char *parName) const
3254 {
3255  if (!parName)
3256  return -1;
3257 
3258  for (Int_t i=0; i<fNpar; i++) {
3259  if (!strcmp(GetParName(i),parName)) return i;
3260  }
3261  return -1;
3262 }
3263 
3264 ////////////////////////////////////////////////////////////////////////////////
3265 /// Return true if the expression at the index 'oper' has to be treated as a string
3266 
3267 Bool_t TFormula::IsString(Int_t oper) const
3268 {
3269  return GetAction(oper) == kStringConst;
3270 }
3271 
3272 ////////////////////////////////////////////////////////////////////////////////
3273 /// Dump this formula with its attributes.
3274 
3275 void TFormula::Print(Option_t *) const
3276 {
3277  Int_t i;
3278  Printf(" %20s : %s Ndim= %d, Npar= %d, Noper= %d",GetName(),GetTitle(), fNdim,fNpar,fNoper);
3279  for (i=0;i<fNoper;i++) {
3280  Printf(" fExpr[%d] = %s action = %d action param = %d ",
3281  i,(const char*)fExpr[i],GetAction(i),GetActionParam(i));
3282  }
3283  //MI change
3284  //
3285  if (fNOperOptimized>0){
3286  Printf("Optimized expression");
3287  for (i=0;i<fNOperOptimized;i++) {
3288  Printf(" fExpr[%d] = %s\t\t action = %d action param = %d ",
3289  i,(const char*)fExprOptimized[i],GetActionOptimized(i),GetActionParamOptimized(i));
3290  }
3291  }
3292 
3293  if (!fNames) return;
3294  if (!fParams) return;
3295  for (i=0;i<fNpar;i++) {
3296  Printf(" Par%3d %20s = %g",i,GetParName(i),fParams[i]);
3297  }
3298 }
3299 
3300 ////////////////////////////////////////////////////////////////////////////////
3301 /// If the formula is for linear fitting, change the title to
3302 /// normal and fill the LinearParts array
3303 
3304 void TFormula::ProcessLinear(TString &formula)
3305 {
3306  TString formula2(formula);
3307  char repl[20];
3308  char *pch;
3309  Int_t nf, offset, replsize;
3310  //replace "++" with "+[i]*"
3311  pch= (char*)strstr(formula.Data(), "++");
3312  if (pch)
3313  formula.Insert(0, "[0]*(");
3314  pch= (char*)strstr(formula.Data(), "++");
3315  if (pch){
3316  //if there are "++", replaces them with +[i]*
3317  nf = 1;
3318  while (pch){
3319  snprintf(repl,20, ")+[%d]*(", nf);
3320  offset = pch-formula.Data();
3321  if (nf<10) replsize = 7;
3322  else if (nf<100) replsize = 8;
3323  else replsize = 9;
3324  formula.Replace(pch-formula.Data(), 2, repl, replsize);
3325  pch = (char*)strstr(formula.Data()+offset, "++");
3326  nf++;
3327  }
3328  formula.Append(')', 1);
3329  } else {
3330  //if there are no ++, create a new string with ++ instead of +[i]*
3331  formula2=formula2(4, formula2.Length()-4);
3332  pch= (char*)strchr(formula2.Data(), '[');
3333  snprintf(repl,20, "++");
3334  nf = 1;
3335  while (pch){
3336  offset = pch-formula2.Data()-1;
3337  if (nf<10) replsize = 5;
3338  else replsize = 6;
3339  formula2.Replace(pch-formula2.Data()-1, replsize, repl, 2);
3340  pch = (char*)strchr(formula2.Data()+offset, '[');
3341  nf++;
3342  }
3343  }
3344 
3345  fLinearParts.Expand(nf);
3346  //break up the formula and fill the array of linear parts
3347  TString replaceformula;
3348  formula2 = formula2.ReplaceAll("++", 2, "|", 1);
3349  TObjArray *oa = formula2.Tokenize("|");
3350  TString replaceformula_name;
3351  for (Int_t i=0; i<nf; i++) {
3352  replaceformula = ((TObjString *)oa->UncheckedAt(i))->GetString();
3353  replaceformula_name = "f_linear_";
3354  replaceformula_name.Append(replaceformula);
3355  TFormula *f = new TFormula(replaceformula_name.Data(), replaceformula.Data());
3356  if (!f) {
3357  Error("TFormula", "f_linear not allocated");
3358  return;
3359  }
3360  {
3361  R__LOCKGUARD(gROOTMutex);
3362  gROOT->GetListOfFunctions()->Remove(f);
3363  }
3364  f->SetBit(kNotGlobal, 1);
3365  fLinearParts.Add(f);
3366  }
3367  oa->Delete();
3368 }
3369 
3370 ////////////////////////////////////////////////////////////////////////////////
3371 /// Initialize parameter number ipar.
3372 
3373 void TFormula::SetParameter(const char *name, Double_t value)
3374 {
3375  Int_t ipar = GetParNumber(name);
3376  if (ipar <0 || ipar >= fNpar) return;
3377  fParams[ipar] = value;
3378  Update();
3379 }
3380 
3381 ////////////////////////////////////////////////////////////////////////////////
3382 /// Initialize parameter number ipar.
3383 
3384 void TFormula::SetParameter(Int_t ipar, Double_t value)
3385 {
3386  if (ipar <0 || ipar >= fNpar) return;
3387  fParams[ipar] = value;
3388  Update();
3389 }
3390 
3391 ////////////////////////////////////////////////////////////////////////////////
3392 /// Initialize array of all parameters.
3393 /// See also the next function with the same name.
3394 
3395 void TFormula::SetParameters(const Double_t *params)
3396 {
3397  for (Int_t i=0; i<fNpar;i++) {
3398  fParams[i] = params[i];
3399  }
3400  Update();
3401 }
3402 
3403 ////////////////////////////////////////////////////////////////////////////////
3404 /// Initialize up to 11 parameters
3405 /// All arguments except THE FIRST TWO are optional
3406 /// In case of a function with only one parameter, call this function with p1=0.
3407 /// Minimum two arguments are required to differentiate this function
3408 /// from the SetParameters(cont Double_t *params)
3409 
3410 void TFormula::SetParameters(Double_t p0,Double_t p1,Double_t p2,Double_t p3,Double_t p4
3411  ,Double_t p5,Double_t p6,Double_t p7,Double_t p8,Double_t p9,Double_t p10)
3412 {
3413  if (fNpar > 0) fParams[0] = p0;
3414  if (fNpar > 1) fParams[1] = p1;
3415  if (fNpar > 2) fParams[2] = p2;
3416  if (fNpar > 3) fParams[3] = p3;
3417  if (fNpar > 4) fParams[4] = p4;
3418  if (fNpar > 5) fParams[5] = p5;
3419  if (fNpar > 6) fParams[6] = p6;
3420  if (fNpar > 7) fParams[7] = p7;
3421  if (fNpar > 8) fParams[8] = p8;
3422  if (fNpar > 9) fParams[9] = p9;
3423  if (fNpar >10) fParams[10]= p10;
3424  Update();
3425 }
3426 
3427 ////////////////////////////////////////////////////////////////////////////////
3428 /// Set name of parameter number ipar
3429 
3430 void TFormula::SetParName(Int_t ipar, const char *name)
3431 {
3432  if (ipar <0 || ipar >= fNpar) return;
3433  fNames[ipar] = name;
3434 }
3435 
3436 ////////////////////////////////////////////////////////////////////////////////
3437 /// Set up to 11 parameter names.
3438 
3439 void TFormula::SetParNames(const char*name0,const char*name1,const char*name2,const char*name3,const char*name4,
3440  const char*name5,const char*name6,const char*name7,const char*name8,const char*name9,const char*name10)
3441 {
3442  if (fNpar > 0) fNames[0] = name0;
3443  if (fNpar > 1) fNames[1] = name1;
3444  if (fNpar > 2) fNames[2] = name2;
3445  if (fNpar > 3) fNames[3] = name3;
3446  if (fNpar > 4) fNames[4] = name4;
3447  if (fNpar > 5) fNames[5] = name5;
3448  if (fNpar > 6) fNames[6] = name6;
3449  if (fNpar > 7) fNames[7] = name7;
3450  if (fNpar > 8) fNames[8] = name8;
3451  if (fNpar > 9) fNames[9] = name9;
3452  if (fNpar >10) fNames[10]= name10;
3453 }
3454 
3455 ////////////////////////////////////////////////////////////////////////////////
3456 /// Stream a class object.
3457 
3458 void TFormula::Streamer(TBuffer &b, const TClass *onfile_class)
3459 {
3460  if (b.IsReading()) {
3461  UInt_t R__s, R__c;
3462  Version_t v = b.ReadVersion(&R__s, &R__c);
3463  if (v==6) {
3464  Error("Streamer","version 6 is not supported");
3465  return;
3466  }
3467  Streamer(b, v, R__s, R__c, onfile_class);
3468 
3469  } else {
3470  b.WriteClassBuffer(TFormula::Class(),this);
3471  }
3472 }
3473 
3474 ////////////////////////////////////////////////////////////////////////////////
3475 /// Stream a class object.
3476 
3477 void TFormula::Streamer(TBuffer &b)
3478 {
3479  if (b.IsReading()) {
3480  UInt_t R__s, R__c;
3481  Version_t v = b.ReadVersion(&R__s, &R__c);
3482  if (v==6) {
3483  Error("Streamer","version 6 is not supported");
3484  return;
3485  }
3486  Streamer(b, v, R__s, R__c, nullptr);
3487 
3488  } else {
3489  b.WriteClassBuffer(TFormula::Class(),this);
3490  }
3491 }
3492 
3493 ////////////////////////////////////////////////////////////////////////////////
3494 /// specialized streamer function being able to read old TF1 versions as TF1Old in memory
3495 
3496 void TFormula::Streamer(TBuffer &b, Int_t v, UInt_t R__s, UInt_t R__c, const TClass *onfile_class)
3497 {
3498 
3499  //printf("Reading TFormula - version %d \n",v);
3500  if (v > 3 ) {
3501  b.ReadClassBuffer(TFormula::Class(), this, v, R__s, R__c, onfile_class);
3502  if (!TestBit(kNotGlobal)) {
3503  R__LOCKGUARD(gROOTMutex);
3504  gROOT->GetListOfFunctions()->Add(this);
3505  }
3506 
3507  // We need to reinstate (if possible) the TMethodCall.
3508  if (fFunctions.GetLast()>=0) {
3509  // Compiles will reset the parameter values so we need
3510  // to temporarily keep them
3511  Double_t *param = fParams;
3512  TString *names = fNames;
3513  Int_t npar = fNpar;
3514  fParams = 0;
3515  fNames = 0;
3516  if (Compile()) {
3517  Error("Streamer","error compiling formula");
3518  return;
3519  }
3520  for (Int_t i = 0; i<npar && i<fNpar; ++i) fParams[i] = param[i];
3521  delete [] param;
3522  delete [] fNames;
3523  fNames = names;
3524  } else if (v<6) {
3525  Convert(v);
3526  }
3527  Optimize();
3528  return;
3529  }
3530  // version smaller or equal to 3
3531  // process old versions before automatic schema evolution
3532  TNamed::Streamer(b);
3533  b >> fNdim;
3534  b >> fNumber;
3535  if (v > 1) b >> fNval;
3536  if (v > 2) b >> fNstring;
3537  fNpar = b.ReadArray(fParams);
3538  fOper = new Int_t[gMAXOP];
3539  fNoper = b.ReadArray(fOper);
3540  fNconst = b.ReadArray(fConst);
3541  if (fNoper) {
3542  fExpr = new TString[fNoper];
3543  }
3544  if (fNpar) {
3545  fNames = new TString[fNpar];
3546  }
3547  Int_t i;
3548  for (i=0;i<fNoper;i++) fExpr[i].Streamer(b);
3549  for (i=0;i<fNpar;i++) fNames[i].Streamer(b);
3550  {
3551  R__LOCKGUARD(gROOTMutex);
3552  if (gROOT->GetListOfFunctions()->FindObject(GetName())) return;
3553  gROOT->GetListOfFunctions()->Add(this);
3554  }
3555  b.CheckByteCount(R__s, R__c, TFormula::IsA());
3556 
3557  Convert(v);
3558  // end of old versions
3559 
3560 }
3561 
3562 void TFormula::Convert(UInt_t /* fromVersion */)
3563 {
3564  // Convert the fOper of a TFormula version fromVersion to the current in memory version
3565 
3566  enum {
3567  kOldexpo = 1000,
3568  kOldgaus = 2000,
3569  kOldlandau = 4000,
3570  kOldxylandau = 4500,
3571  kOldConstants = 50000,
3572  kOldStrings = 80000,
3573  kOldVariable = 100000,
3574  kOldTreeString = 105000,
3575  kOldFormulaVar = 110000,
3576  kOldBoolOptimize = 120000,
3577  kOldFunctionCall = 200000
3578  };
3579  int i,j;
3580 
3581  for (i=0,j=0; i<fNoper; ++i,++j) {
3582  Int_t action = fOper[i];
3583  Int_t newActionCode = 0;
3584  Int_t newActionParam = 0;
3585 
3586  if ( action == 0) {
3587  // Sign Inversion
3588 
3589  newActionCode = kSignInv;
3590 
3591  Float_t aresult = 99.99;
3592  sscanf((const char*)fExpr[i],"%g",&aresult);
3593  R__ASSERT((aresult+1)<0.001);
3594 
3595  ++i; // skip the implied multiplication.
3596 
3597  // For consistency and for Optimize to work correctly
3598  // we need to remove the "-1" string in fExpr
3599  for (int z=i; z<fNoper; ++z) {
3600  fExpr[z-1] = fExpr[z];
3601  }
3602 
3603  } else if ( action < 100 ) {
3604  // basic operators and mathematical library
3605 
3606  newActionCode = action;
3607 
3608  } else if (action >= kOldFunctionCall) {
3609  // Function call
3610 
3611  newActionCode = kFunctionCall;
3612  newActionParam = action-kOldFunctionCall;
3613 
3614  } else if (action >= kOldBoolOptimize) {
3615  // boolean operation optimizer
3616 
3617  newActionCode = kBoolOptimize;
3618  newActionParam = action-kOldBoolOptimize;
3619 
3620  } else if (action >= kOldFormulaVar) {
3621  // a variable
3622 
3623  newActionCode = kVariable;
3624  newActionParam = action-kOldFormulaVar;
3625 
3626  } else if (action >= kOldTreeString) {
3627  // a tree string
3628 
3629  newActionCode = kDefinedString;
3630  newActionParam = action-kOldTreeString;
3631 
3632  } else if (action >= kOldVariable) {
3633  // a tree variable
3634 
3635  newActionCode = kDefinedVariable;
3636  newActionParam = action-kOldVariable;
3637 
3638  } else if (action == kOldStrings) {
3639  // String
3640 
3641  newActionCode = kStringConst;
3642 
3643  } else if (action >= kOldConstants) {
3644  // numerical value
3645 
3646  newActionCode = kConstant;
3647  newActionParam = action-kOldConstants;
3648 
3649  } else if (action > 10000 && action < kOldConstants) {
3650  // Polynomial
3651 
3652  int var = action/10000; //arrondit
3653  newActionCode = kpol + (var-1);
3654  newActionParam = action - var*10000;
3655 
3656  } else if (action >= 4600) {
3657 
3658  Error("Convert","Unsupported value %d",action);
3659 
3660  } else if (action > kOldxylandau) {
3661  // xylandau
3662 
3663  newActionCode = kxylandau;
3664  newActionParam = action - (kOldxylandau+1);
3665 
3666  } else if (action > kOldlandau) {
3667  // landau, xlandau, ylandau or zlandau
3668 
3669  newActionCode = klandau;
3670  int var = action/100-40;
3671  if (var) newActionCode += var;
3672  newActionParam = action - var*100 - (kOldlandau+1);
3673 
3674  } else if (action > 2500 && action < 2600) {
3675  // xygaus
3676 
3677  newActionCode = kxygaus;
3678  newActionParam = action-2501;
3679 
3680  } else if (action > 2000 && action < 2500) {
3681  // gaus, xgaus, ygaus or zgaus
3682 
3683  newActionCode = kgaus;
3684  int var = action/100-20;
3685  if (var) newActionCode += var;
3686  newActionParam = action - var*100 - (kOldgaus+1);
3687 
3688  } else if (action > 1500 && action < 1600) {
3689  // xyexpo
3690 
3691  newActionCode = kxyexpo;
3692  newActionParam = action-1501;
3693 
3694  } else if (action > 1000 && action < 1500) {
3695  // expo or xexpo or yexpo or zexpo
3696 
3697  newActionCode = kexpo;
3698  int var = action/100-10;
3699  if (var) newActionCode += var;
3700  newActionParam = action - var*100 - (kOldexpo+1);
3701 
3702  } else if (action > 100 && action < 200) {
3703  // Parameter substitution
3704 
3705  newActionCode = kParameter;
3706  newActionParam = action - 101;
3707  }
3708 
3709  SetAction( j, newActionCode, newActionParam );
3710 
3711  }
3712  if (i!=j) {
3713  fNoper -= (i-j);
3714  }
3715 
3716 }
3717 
3718 ////////////////////////////////////////////////////////////////////////////////
3719 /// TOper offset - helper class for TFormula*
3720 /// specify type of operand
3721 /// fTypeX = kVariable
3722 /// = kParameter
3723 /// = kConstant
3724 /// fOffestX = offset in corresponding array
3725 
3726 TOperOffset::TOperOffset()
3727 {
3728  fType0=0;
3729  fType1=0;
3730  fType2=0;
3731  fType3=0;
3732  fOffset0=0;
3733  fOffset1=0;
3734  fOffset2=0;
3735  fOffset3=0;
3736  fOldAction=0;
3737  fToJump=0;
3738 }
3739 
3740 ////////////////////////////////////////////////////////////////////////////////
3741 /// MakePrimitive
3742 /// find TFormulaPrimitive replacement for some operands
3743 
3744 void TFormula::MakePrimitive(const char *expr, Int_t pos)
3745 {
3746  TString cbase(expr);
3747  cbase.ReplaceAll("Double_t ","");
3748  int paran = cbase.First("(");
3749  // int nargs = 0;
3750  if (paran>0) {
3751  //nargs = 1;
3752  cbase[paran]=0;
3753  }
3754 
3755  if (cbase=="<") cbase="XlY";
3756  if (cbase=="<=") cbase="XleY";
3757  if (cbase==">") cbase="XgY";
3758  if (cbase==">=") cbase="XgeY";
3759  if (cbase=="==" && GetActionOptimized(pos)!=kStringEqual) cbase="XeY";
3760  if (cbase=="!=" && GetActionOptimized(pos)!=kStringNotEqual) cbase="XneY";
3761 
3762  ROOT::v5::TFormulaPrimitive *prim = ROOT::v5::TFormulaPrimitive::FindFormula(cbase ,paran>0 ? cbase.Data() + paran + 1 : (const char*)0);
3763  if (prim) {
3764  fPredefined[pos] = prim;
3765  if (prim->fType==10) {
3766  SetActionOptimized(pos, kFD1);
3767  }
3768  if (prim->fType==110) {
3769  SetActionOptimized(pos, kFD2);
3770  }
3771  if (prim->fType==1110) {
3772  SetActionOptimized(pos, kFD3);
3773  }
3774  if (prim->fType==-1) {
3775  SetActionOptimized(pos, kFDM);
3776  }
3777  if (prim->fType==0){
3778  SetActionOptimized(pos,kConstant,fNconst);
3779  fConst[fNconst] = prim->Eval(0);
3780  fNconst++;
3781  }
3782  return;
3783  }
3784 }
3785 
3786 ////////////////////////////////////////////////////////////////////////////////
3787 /// MI include
3788 ///
3789 /// Optimize formula
3790 /// - Minimize the number of operands
3791 /// 1. several operands are glued together
3792 /// 2. some primitive functions glued together - exemp. (x+y) => PlusXY(x,y)
3793 /// 3. maximize number of standard calls minimizing number of jumps in Eval cases
3794 /// 4. variables, parameters and constants are mapped - using fOperOfssets0
3795 /// Eval procedure use direct acces to data (only one corresponding case statement in eval procedure)
3796 /// ~~~ {.cpp}
3797 /// pdata[operand={Var,Par,Const}][offset]
3798 /// pdata[fOperOffsets0[i]][fOperOffset1[i+1]]
3799 /// ~~~
3800 /// - The fastest evaluation function is chosen at the end
3801 /// 1. fOptimal := pointer to the fastest function for given evaluation string
3802 /// ~~~ {.cpp}
3803 /// switch(GetActionOptimized(0)){
3804 /// case kData : {fOptimal= (TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive0; break;}
3805 /// case kUnary : {fOptimal= (TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive1; break;}
3806 /// case kBinary : {fOptimal= (TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive2; break;}
3807 /// case kThree : {fOptimal= (TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive3; break;}
3808 /// case kFDM : {fOptimal= (TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive4; break;}
3809 /// }
3810 /// ~~~
3811 /// 2. ex.
3812 /// - fOptimal = ::EvalPrimitive0 - if it return only variable, constant or parameter
3813 /// - = ::EvalParameter1 - if only one unary operation
3814 /// - = ::EvalPrimitive2 - if only one binary operation
3815 
3816 void TFormula::Optimize()
3817 {
3818  //
3819  // Initialize data members
3820  //
3821 
3822  Int_t i;
3823 
3824  if (fPredefined) { delete [] fPredefined; fPredefined = 0;}
3825  if (fOperOffset) { delete [] fOperOffset; fOperOffset = 0;}
3826  if (fExprOptimized) { delete [] fExprOptimized; fExprOptimized = 0;}
3827  if (fOperOptimized) { delete [] fOperOptimized; fOperOptimized = 0;}
3828 
3829  fExprOptimized = new TString[fNoper];
3830  fOperOptimized = new Int_t[fNoper];
3831  fPredefined = new ROOT::v5::TFormulaPrimitive*[fNoper];
3832  fOperOffset = new TOperOffset[fNoper];
3833  for (i=0; i<fNoper; i++) {
3834  fExprOptimized[i] = fExpr[i] ;
3835  fOperOptimized[i] = fOper[i];
3836  fPredefined[i]= 0;
3837  }
3838 
3839  //
3840  //Make primitives
3841  //
3842  for (i=0;i<fNoper;i++){
3843  if (fExprOptimized[i].Data()) {
3844  MakePrimitive(fExprOptimized[i].Data(), i);
3845  }
3846  }
3847  //
3848  Int_t maxfound = fNoper+1;
3849  Int_t *offset = new Int_t[maxfound*16];
3850  Int_t *optimized = new Int_t[maxfound];
3851  //
3852  //
3853  ROOT::v5::TFormulaPrimitive* primitive[10];
3854  primitive[0] = ROOT::v5::TFormulaPrimitive::FindFormula("PlusXY");
3855  primitive[1] = ROOT::v5::TFormulaPrimitive::FindFormula("MinusXY");
3856  primitive[2] = ROOT::v5::TFormulaPrimitive::FindFormula("MultXY");
3857  primitive[3] = ROOT::v5::TFormulaPrimitive::FindFormula("DivXY");
3858  primitive[4] = ROOT::v5::TFormulaPrimitive::FindFormula("XpYpZ");
3859  primitive[5] = ROOT::v5::TFormulaPrimitive::FindFormula("XxYxZ");
3860  primitive[6] = ROOT::v5::TFormulaPrimitive::FindFormula("XxYpZ");
3861  primitive[7] = ROOT::v5::TFormulaPrimitive::FindFormula("XpYxZ");
3862  primitive[8] = ROOT::v5::TFormulaPrimitive::FindFormula("Pow2");
3863  primitive[9] = ROOT::v5::TFormulaPrimitive::FindFormula("Pow3");
3864  //
3865  // set data pointers
3866  //
3867  for (i=0;i<fNoper;i++) optimized[i]=0;
3868  //
3869  for (i=0;i<fNoper;i++){
3870  Int_t actionparam = GetActionParamOptimized(i);
3871  Int_t action = GetActionOptimized(i);
3872 
3873  if (action==kBoolOptimize){
3874  //
3875  // optimize booleans
3876  //
3877  fOperOffset[i].fType1 = actionparam/10; // operands to skip
3878  fOperOffset[i].fOffset0 = actionparam%10; // 1 is && , 2 is || - operand
3879  fOperOffset[i].fToJump = i+fOperOffset[i].fType1; // where we should jump
3880  continue;
3881  }
3882  if (action==kJump || action==kJumpIf) {
3883  // Ternary conditional operator
3884  fOperOffset[i].fType1 = action;
3885  fOperOffset[i].fToJump = actionparam;
3886  }
3887  //
3888  if (action==kConstant&&i<fNoper-2){
3889  //
3890  // get offsets for kFDM operands
3891  //
3892  if (GetActionOptimized(i+1)==kConstant && GetActionOptimized(i+2)==kFDM){
3893  optimized[i]=1;
3894  optimized[i+1]=1;
3895  i+=2;
3896  fOperOffset[i].fType0=actionparam;
3897  fOperOffset[i].fOffset0=GetActionParamOptimized(i-1);
3898  Int_t offset2 = int(fConst[fOperOffset[i].fOffset0]+0.4);
3899  fOperOffset[i].fOffset0=offset2;
3900  Int_t nparmax = offset2+fPredefined[i]->fNParameters;
3901  if (nparmax>fNpar){ // increase expected number of parameters
3902  fNpar=nparmax;
3903  }
3904  continue;
3905  }
3906  }
3907  switch(action){
3908  case kVariable : {action=kData; fOperOffset[i].fType0=0; break;}
3909  case kParameter: {action=kData; fOperOffset[i].fType0=1; break;}
3910  case kConstant : {action=kData; fOperOffset[i].fType0=2; break;}
3911  }
3912  //
3913  fOperOffset[i].fOffset0 = GetActionParamOptimized(i);
3914  SetActionOptimized(i,action, actionparam); //set common data option
3915  }
3916  //
3917  //
3918  fNOperOptimized = fNoper;
3919  //
3920  for (i=0; i<fNoper; ++i)
3921  {
3922  //
3923  if (!(GetActionOptimized(i)== kData)) continue;
3924  offset[0] = fOperOffset[i].fType0; //
3925  offset[1] = fOperOptimized[i] & kTFOperMask; // offset
3926 
3927  if ((i+1) >= fNoper) continue;
3928 
3929  if (GetActionOptimized(i+1)==kFD1){
3930  optimized[i] = 1; // to be optimized
3931  i++;
3932  fOperOffset[i].fType0 = offset[0];
3933  fOperOffset[i].fOffset0 = offset[1];
3934  SetActionOptimized(i ,kUnary);
3935  continue;
3936  }
3937  if (GetActionOptimized(i+1)==kAdd){
3938  optimized[i] = 1; // to be optimized
3939  i++;
3940  fOperOffset[i].fType0 = offset[0];
3941  fOperOffset[i].fOffset0 = offset[1];
3942  SetActionOptimized(i ,kPlusD);
3943  continue;
3944  }
3945  if (GetActionOptimized(i+1)==kMultiply){
3946  optimized[i] = 1; // to be optimized
3947  i++;
3948  fOperOffset[i].fType0 = offset[0];
3949  fOperOffset[i].fOffset0 = offset[1];
3950  SetActionOptimized(i,kMultD);
3951  continue;
3952  }
3953 
3954  if ((i+2) >= fNoper) continue;
3955 
3956  //
3957  //Binary operators
3958  if (!(GetActionOptimized(i+1)== kData)) continue;
3959  offset[2] = fOperOffset[i+1].fType0;
3960  offset[3] = fOperOptimized[i+1] & kTFOperMask; // offset
3961  //
3962  if (GetActionOptimized(i+2)==kFD2 || GetActionOptimized(i+2)==kAdd ||GetActionOptimized(i+2)==kSubstract||
3963  GetActionOptimized(i+2)==kMultiply || GetActionOptimized(i+2)==kDivide){
3964 
3965  optimized[i] = 1; // to be optimized
3966  optimized[i+1] = 1; // to be optimized
3967  i+=2;
3968  //
3969  fOperOffset[i].fType0 = offset[0];
3970  fOperOffset[i].fOffset0 = offset[1];
3971  fOperOffset[i].fType1 = offset[2];
3972  fOperOffset[i].fOffset1 = offset[3];
3973  fOperOffset[i].fType2 = GetActionOptimized(i); //remember old action
3974  if (GetActionOptimized(i)==kAdd) {fPredefined[i] = primitive[0];}
3975  if (GetActionOptimized(i)==kSubstract) {fPredefined[i] = primitive[1];}
3976  if (GetActionOptimized(i)==kMultiply) {
3977  fPredefined[i]=primitive[2];
3978  if (offset[0]==offset[2]&&offset[1]==offset[3]) {
3979  fPredefined[i] = primitive[8];
3980  SetActionOptimized(i,kUnary);
3981  continue;
3982  }
3983  }
3984  if (GetActionOptimized(i)==kDivide) {
3985  fPredefined[i] = primitive[3];
3986  }
3987  SetActionOptimized(i,kBinary);
3988  continue;
3989  }
3990 
3991  if ((i+3) >= fNoper) continue;
3992 
3993  //
3994  //operator 3
3995  //
3996  if (!(GetActionOptimized(i+2)== kData)) continue;
3997  offset[4] = fOperOffset[i+2].fType0;
3998  offset[5] = fOperOptimized[i+2] & kTFOperMask; // offset
3999  //
4000  if (GetActionOptimized(i+3)==kFD3|| ( (GetActionOptimized(i+3)==kAdd||GetActionOptimized(i+3)==kMultiply) &&
4001  (GetActionOptimized(i+4)==kAdd||GetActionOptimized(i+4)==kMultiply) ) ){
4002  optimized[i+0] = 1; // to be optimized
4003  optimized[i+1] = 1; // to be optimized
4004  optimized[i+2] = 1; // to be optimized
4005  i+=3;
4006  //
4007  fOperOffset[i].fType0 = offset[0];
4008  fOperOffset[i].fOffset0 = offset[1];
4009  fOperOffset[i].fType1 = offset[2];
4010  fOperOffset[i].fOffset1 = offset[3];
4011  fOperOffset[i].fType2 = offset[4];
4012  fOperOffset[i].fOffset2 = offset[5];
4013  //
4014  fOperOffset[i].fOldAction = GetActionOptimized(i); //remember old action
4015  if (GetActionOptimized(i)==kFD3) {
4016  SetActionOptimized(i,kThree);
4017  continue;
4018  }
4019  Int_t action=0;
4020  Int_t action2=kThree;
4021  if (GetActionOptimized(i)==kAdd&&GetActionOptimized(i+1)==kAdd) action=4;
4022  if (GetActionOptimized(i)==kMultiply&&GetActionOptimized(i+1)==kMultiply) {
4023  action=5;
4024  if (offset[0]==offset[2]&&offset[1]==offset[3]&&offset[0]==offset[4]&&offset[1]==offset[5]){
4025  fPredefined[i]=primitive[9];
4026  action2=kUnary;
4027  action =9;
4028  }
4029  }
4030  if (GetActionOptimized(i)==kAdd&&GetActionOptimized(i+1)==kMultiply) action=6;
4031  if (GetActionOptimized(i)==kMultiply&&GetActionOptimized(i+1)==kAdd) action=7;
4032  //
4033  optimized[i]=1;
4034  i++;
4035  fOperOffset[i].fType0 = offset[0];
4036  fOperOffset[i].fOffset0 = offset[1];
4037  fOperOffset[i].fType1 = offset[2];
4038  fOperOffset[i].fOffset1 = offset[3];
4039  fOperOffset[i].fType2 = offset[4];
4040  fOperOffset[i].fOffset2 = offset[5];
4041  fPredefined[i]=primitive[action];
4042  SetActionOptimized(i,action2);
4043  continue;
4044  }
4045  }
4046  //
4047  //
4048  Int_t operO=0;
4049  TString expr="";
4050  Int_t *map0 = new Int_t[maxfound]; //remapping of the operands
4051  Int_t *map1 = new Int_t[maxfound]; //remapping of the operands
4052  for (i=0;i<fNoper;i++){
4053  map0[i] = operO;
4054  map1[operO] = i;
4055  fOperOptimized[operO] = fOperOptimized[i];
4056  fPredefined[operO] = fPredefined[i];
4057  fOperOffset[operO] = fOperOffset[i];
4058  expr += fExprOptimized[i];
4059  if (optimized[i]==0){
4060  fExprOptimized[operO] = expr;
4061  expr = "";
4062  operO++;
4063  }else{
4064  expr += ",";
4065  }
4066  }
4067  //
4068  // Recalculate long jump for Boolean optimize
4069  //
4070  for (i=0; i<fNOperOptimized; i++){
4071  Int_t optaction = GetActionOptimized(i);
4072  if (optaction==kBoolOptimize){
4073  Int_t oldpos = fOperOffset[i].fToJump;
4074  Int_t newpos = oldpos==fNoper ? fNOperOptimized : map0[oldpos];
4075  fOperOffset[i].fToJump = newpos; // new position to jump
4076  Int_t actionop = GetActionParamOptimized(i) % 10;
4077  switch (actionop) {
4078  case 1: SetActionOptimized(i,kBoolOptimizeAnd,newpos); break;
4079  case 2: SetActionOptimized(i,kBoolOptimizeOr,newpos); break;
4080  }
4081  } else if (optaction==kJump || optaction==kJumpIf) {
4082  Int_t oldpos = fOperOffset[i].fToJump;
4083  Int_t newpos = oldpos==fNoper ? fNOperOptimized : map0[oldpos];
4084  fOperOffset[i].fToJump = newpos; // new position to jump
4085  SetActionOptimized(i,optaction,newpos);
4086  }
4087  }
4088 
4089 
4090  fNOperOptimized = operO;
4091  //
4092  fOptimal= (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalParFast;
4093  if (fNOperOptimized==1) {
4094  switch(GetActionOptimized(0)){
4095  case kData : {fOptimal= (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive0; break;}
4096  case kUnary : {fOptimal= (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive1; break;}
4097  case kBinary : {fOptimal= (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive2; break;}
4098  case kThree : {fOptimal= (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive3; break;}
4099  case kFDM : {fOptimal= (ROOT::v5::TFormulaPrimitive::TFuncG)&TFormula::EvalPrimitive4; break;}
4100  }
4101  }
4102 
4103  delete [] map1;
4104  delete [] map0;
4105  delete [] offset;
4106  delete [] optimized;
4107 }
4108 
4109 ////////////////////////////////////////////////////////////////////////////////
4110 /// Evaluate primitive formula
4111 
4112 Double_t TFormula::EvalPrimitive(const Double_t *x, const Double_t *params)
4113 {
4114  const Double_t *pdata[3] = {x,(params!=0)?params:fParams, fConst};
4115  Double_t result = pdata[fOperOffset->fType0][fOperOffset->fOffset0];
4116  switch((fOperOptimized[0] >> kTFOperShift)) {
4117  case kData : return result;
4118  case kUnary : return (fPredefined[0]->fFunc10)(pdata[fOperOffset->fType0][fOperOffset->fOffset0]);
4119  case kBinary :return (fPredefined[0]->fFunc110)(result,
4120  pdata[fOperOffset->fType1][fOperOffset->fOffset1]);
4121 
4122  case kThree :return (fPredefined[0]->fFunc1110)(result, pdata[fOperOffset->fType1][fOperOffset->fOffset1],
4123  pdata[fOperOffset->fType2][fOperOffset->fOffset2]);
4124  case kFDM : return (fPredefined[0]->fFuncG)((Double_t*)&x[fOperOffset->fType0],
4125  (Double_t*)&params[fOperOffset->fOffset0]);
4126  }
4127  return 0;
4128 }
4129 
4130 ////////////////////////////////////////////////////////////////////////////////
4131 /// Evaluate primitive formula
4132 
4133 Double_t TFormula::EvalPrimitive0(const Double_t *x, const Double_t *params)
4134 {
4135  const Double_t *pdata[3] = {x,(params!=0)?params:fParams, fConst};
4136  return pdata[fOperOffset->fType0][fOperOffset->fOffset0];
4137 }
4138 
4139 ////////////////////////////////////////////////////////////////////////////////
4140 /// Evaluate primitive formula
4141 
4142 Double_t TFormula::EvalPrimitive1(const Double_t *x, const Double_t *params)
4143 {
4144  const Double_t *pdata[3] = {x,(params!=0)?params:fParams, fConst};
4145  return (fPredefined[0]->fFunc10)(pdata[fOperOffset->fType0][fOperOffset->fOffset0]);
4146 }
4147 
4148 ////////////////////////////////////////////////////////////////////////////////
4149 /// Evaluate primitive formula
4150 
4151 Double_t TFormula::EvalPrimitive2(const Double_t *x, const Double_t *params)
4152 {
4153  const Double_t *pdata[3] = {x,(params!=0)?params:fParams, fConst};
4154  return (fPredefined[0]->fFunc110)(pdata[fOperOffset->fType0][fOperOffset->fOffset0],
4155  pdata[fOperOffset->fType1][fOperOffset->fOffset1]);
4156 }
4157 
4158 ////////////////////////////////////////////////////////////////////////////////
4159 /// Evaluate primitive formula
4160 
4161 Double_t TFormula::EvalPrimitive3(const Double_t *x, const Double_t *params)
4162 {
4163  const Double_t *pdata[3] = {x,(params!=0)?params:fParams, fConst};
4164  return (fPredefined[0]->fFunc1110)(pdata[fOperOffset->fType0][fOperOffset->fOffset0], pdata[fOperOffset->fType1][fOperOffset->fOffset1],
4165  pdata[fOperOffset->fType2][fOperOffset->fOffset2]);
4166 }
4167 
4168 ////////////////////////////////////////////////////////////////////////////////
4169 /// Evaluate primitive formula
4170 
4171 Double_t TFormula::EvalPrimitive4(const Double_t *x, const Double_t *params)
4172 {
4173  const Double_t *par = (params!=0)?params:fParams;
4174  return (fPredefined[0]->fFuncG)((Double_t*)&x[fOperOffset->fType0],
4175  (Double_t*)&par[fOperOffset->fOffset0]);
4176 }
4177 
4178 ////////////////////////////////////////////////////////////////////////////////
4179 /// Evaluate this formula.
4180 ///
4181 /// The current value of variables x,y,z,t is passed through the pointer x.
4182 /// The parameters used will be the ones in the array params if params is given
4183 /// otherwise parameters will be taken from the stored data members fParams
4184 ///
4185 /// \image html TFormula_eval.png
4186 
4187 Double_t TFormula::EvalParFast(const Double_t *x, const Double_t *uparams)
4188 {
4189  const Double_t *pdata[3] = {x,(uparams!=0)?uparams:fParams, fConst};
4190  //
4191  Int_t i,j;
4192  Double_t tab[kMAXFOUND] = {0};
4193  const char *stringStack[gMAXSTRINGFOUND] = {0};
4194  Double_t param_calc[kMAXFOUND];
4195  char *string_calc[gMAXSTRINGFOUND] = {0};
4196  Int_t precalculated = 0;
4197  Int_t precalculated_str = 0;
4198 
4199  Double_t *params;
4200 
4201  if (uparams) {
4202  //for (j=0;j<fNpar;j++) fParams[j] = params[j];
4203  params = const_cast<Double_t*>(uparams);
4204  } else {
4205  params = fParams;
4206  }
4207 
4208  //if (params) {
4209  // for (j=0;j<fNpar;j++) fParams[j] = params[j];
4210  //}
4211  UInt_t pos = 0;
4212  UInt_t strpos = 0;
4213  // for (i=0; i<fNoper; ++i) {
4214  for (i=0; i<fNOperOptimized; ++i) {
4215  //
4216  const int oper = fOperOptimized[i];
4217  const int opcode = oper >> kTFOperShift;
4218 
4219  switch(opcode) { // FREQUENTLY USED OPERATION
4220  case kData : tab[pos] = pdata[fOperOffset[i].fType0][fOperOffset[i].fOffset0]; pos++;continue;
4221  case kPlusD : tab[pos-1]+= pdata[fOperOffset[i].fType0][fOperOffset[i].fOffset0]; continue;
4222  case kMultD : tab[pos-1]*= pdata[fOperOffset[i].fType0][fOperOffset[i].fOffset0]; continue;
4223  case kAdd : pos--; tab[pos-1] += tab[pos]; continue;
4224  case kSubstract : pos--; tab[pos-1] -= tab[pos]; continue;
4225  case kMultiply : pos--; tab[pos-1] *= tab[pos]; continue;
4226  case kDivide : pos--; if (tab[pos] == 0) tab[pos-1] = 0; // division by 0
4227  else tab[pos-1] /= tab[pos];
4228  continue;
4229  case kUnary : tab[pos] = (fPredefined[i]->fFunc10)(pdata[fOperOffset[i].fType0][fOperOffset[i].fOffset0]); pos++;continue;
4230  case kBinary : tab[pos] = (fPredefined[i]->fFunc110)(pdata[fOperOffset[i].fType0][fOperOffset[i].fOffset0],
4231  pdata[fOperOffset[i].fType1][fOperOffset[i].fOffset1]);pos++;continue;
4232 
4233  case kThree : tab[pos] = (fPredefined[i]->fFunc1110)(pdata[fOperOffset[i].fType0][fOperOffset[i].fOffset0],
4234  pdata[fOperOffset[i].fType1][fOperOffset[i].fOffset1],
4235  pdata[fOperOffset[i].fType2][fOperOffset[i].fOffset2]); pos++; continue;
4236 
4237  case kFDM : tab[pos] = (fPredefined[i]->fFuncG)(&x[fOperOffset[i].fType0],&params[fOperOffset[i].fOffset0]); pos++;continue;
4238  case kFD1 : tab[pos-1] =(fPredefined[i]->fFunc10)(tab[pos-1]); continue;
4239  case kFD2 : pos--; tab[pos-1] = (fPredefined[i]->fFunc110)(tab[pos-1],tab[pos]); continue;
4240  case kFD3 : pos-=2; tab[pos-1] = (fPredefined[i]->fFunc1110)(tab[pos-1],tab[pos],tab[pos+1]); continue;
4241  }
4242  //
4243  switch(opcode) {
4244  case kBoolOptimizeAnd: {
4245  if (!tab[pos-1]) i=fOperOffset[i].fToJump;
4246  continue;
4247  }
4248  case kBoolOptimizeOr: {
4249  if (tab[pos-1]) i=fOperOffset[i].fToJump;
4250  continue;
4251  }
4252  case kAnd : pos--; tab[pos-1] = (bool)tab[pos]; continue; // use the fact that other were check before - see bool optimize
4253  case kOr : pos--; tab[pos-1] = (bool)tab[pos]; continue;
4254  }
4255  switch(opcode) {
4256  // case kabs : tab[pos-1] = TMath::Abs(tab[pos-1]); continue;
4257  case kabs : if (tab[pos-1]<0) tab[pos-1]=-tab[pos-1]; continue;
4258  case ksign : if (tab[pos-1] < 0) tab[pos-1] = -1; else tab[pos-1] = 1; continue;
4259  case kint : tab[pos-1] = Double_t(Int_t(tab[pos-1])); continue;
4260  case kpow : pos--; tab[pos-1] = TMath::Power(tab[pos-1],tab[pos]); continue;
4261 
4262  case kModulo : {pos--;
4263  Long64_t int1((Long64_t)tab[pos-1]);
4264  Long64_t int2((Long64_t)tab[pos]);
4265  tab[pos-1] = Double_t(int1%int2);
4266  continue;}
4267 
4268 
4269  case kStringConst: { strpos++; stringStack[strpos-1] = (char*)fExprOptimized[i].Data(); pos++; tab[pos-1] = 0; continue; }
4270  case kfmod : pos--; tab[pos-1] = fmod(tab[pos-1],tab[pos]); continue;
4271 
4272  case kstrstr : strpos -= 2; pos-=2; pos++;
4273  if (strstr(stringStack[strpos],stringStack[strpos+1])) tab[pos-1]=1;
4274  else tab[pos-1]=0;
4275  continue;
4276  case kpi : pos++; tab[pos-1] = TMath::ACos(-1); continue;
4277 
4278 
4279  case kSignInv: tab[pos-1] = -1 * tab[pos-1]; continue;
4280 
4281  case krndm : pos++; tab[pos-1] = gRandom->Rndm(); continue;
4282 
4283 
4284  case kEqual: pos--; if (tab[pos-1] == tab[pos]) tab[pos-1]=1;
4285  else tab[pos-1]=0;
4286  continue;
4287  case kNotEqual : pos--; if (tab[pos-1] != tab[pos]) tab[pos-1]=1;
4288  else tab[pos-1]=0;
4289  continue;
4290  case kNot : if (tab[pos-1]!=0) tab[pos-1] = 0; else tab[pos-1] = 1;
4291  continue;
4292 
4293  case kStringEqual : strpos -= 2; pos -=2 ; pos++;
4294  if (!strcmp(stringStack[strpos+1],stringStack[strpos])) tab[pos-1]=1;
4295  else tab[pos-1]=0;
4296  continue;
4297  case kStringNotEqual: strpos -= 2; pos -= 2; pos++;
4298  if (strcmp(stringStack[strpos+1],stringStack[strpos])) tab[pos-1]=1;
4299  else tab[pos-1]=0;
4300  continue;
4301 
4302  case kBitAnd : pos--; tab[pos-1]= ((Int_t) tab[pos-1]) & ((Int_t) tab[pos]); continue;
4303  case kBitOr : pos--; tab[pos-1]= ((Int_t) tab[pos-1]) | ((Int_t) tab[pos]); continue;
4304  case kLeftShift : pos--; tab[pos-1]= ((Int_t) tab[pos-1]) <<((Int_t) tab[pos]); continue;
4305  case kRightShift: pos--; tab[pos-1]= ((Int_t) tab[pos-1]) >>((Int_t) tab[pos]); continue;
4306 
4307  case kJump : i = (oper & kTFOperMask); continue;
4308  case kJumpIf : pos--; if (!tab[pos]) i = (oper & kTFOperMask); continue;
4309 
4310  case kBoolOptimize: {
4311  // boolean operation optimizer
4312 
4313  int param = (oper & kTFOperMask);
4314  int op = param % 10; // 1 is && , 2 is ||
4315 
4316  if (op == 1 && (!tab[pos-1]) ) {
4317  // &&: skip the right part if the left part is already false
4318 
4319  i += param / 10;
4320 
4321  // Preserve the existing behavior (i.e. the result of a&&b is
4322  // either 0 or 1)
4323  tab[pos-1] = 0;
4324 
4325  } else if (op == 2 && tab[pos-1] ) {
4326  // ||: skip the right part if the left part is already true
4327 
4328  i += param / 10;
4329 
4330  // Preserve the existing behavior (i.e. the result of a||b is
4331  // either 0 or 1)
4332  tab[pos-1] = 1;
4333 
4334  }
4335 
4336  continue;
4337  }
4338 
4339  }
4340  switch(opcode) {
4341 
4342 #define R__EXPO(var) \
4343  { \
4344  pos++; int param = (oper & kTFOperMask); \
4345  tab[pos-1] = TMath::Exp(params[param]+params[param+1]*x[var]); \
4346  continue; \
4347  }
4348  // case kexpo:
4349  case kxexpo: R__EXPO(0);
4350  case kyexpo: R__EXPO(1);
4351  case kzexpo: R__EXPO(2);
4352  case kxyexpo:{ pos++; int param = (oper & kTFOperMask);
4353  tab[pos-1] = TMath::Exp(params[param]+params[param+1]*x[0]+params[param+2]*x[1]);
4354  continue; }
4355 #ifdef R__GAUS
4356 #undef R__GAUS
4357 #endif
4358 #define R__GAUS(var) \
4359  { \
4360  pos++; int param = (oper & kTFOperMask); \
4361  tab[pos-1] = params[param]*TMath::Gaus(x[var],params[param+1], \
4362  params[param+2],IsNormalized()); \
4363  continue; \
4364  }
4365 
4366  // case kgaus:
4367  case kxgaus: R__GAUS(0);
4368  case kygaus: R__GAUS(1);
4369  case kzgaus: R__GAUS(2);
4370  case kxygaus: { pos++; int param = (oper & kTFOperMask);
4371  Double_t intermede1;
4372  if (params[param+2] == 0) {
4373  intermede1=1e10;
4374  } else {
4375  intermede1=Double_t((x[0]-params[param+1])/params[param+2]);
4376  }
4377  Double_t intermede2;
4378  if (params[param+4] == 0) {
4379  intermede2=1e10;
4380  } else {
4381  intermede2=Double_t((x[1]-params[param+3])/params[param+4]);
4382  }
4383  tab[pos-1] = params[param]*TMath::Exp(-0.5*(intermede1*intermede1+intermede2*intermede2));
4384  continue; }
4385 
4386 #define R__LANDAU(var) \
4387  { \
4388  pos++; const int param = (oper & kTFOperMask); \
4389  tab[pos-1] = params[param]*TMath::Landau(x[var],params[param+1],params[param+2],IsNormalized()); \
4390  continue; \
4391  }
4392  // case klandau:
4393  case kxlandau: R__LANDAU(0);
4394  case kylandau: R__LANDAU(1);
4395  case kzlandau: R__LANDAU(2);
4396  case kxylandau: { pos++; int param = oper&0x7fffff /* ActionParams[i] */ ;
4397  Double_t intermede1=TMath::Landau(x[0], params[param+1], params[param+2],IsNormalized());
4398  Double_t intermede2=TMath::Landau(x[1], params[param+3], params[param+4],IsNormalized());
4399  tab[pos-1] = params[param]*intermede1*intermede2;
4400  continue;
4401  }
4402 
4403 #define R__POLY(var) \
4404  { \
4405  pos++; int param = (oper & kTFOperMask); \
4406  tab[pos-1] = 0; Double_t intermede = 1; \
4407  Int_t inter = param/100; /* arrondit */ \
4408  Int_t int1= param-inter*100-1; /* aucune simplification ! (sic) */ \
4409  for (j=0 ;j<inter+1;j++) { \
4410  tab[pos-1] += intermede*params[j+int1]; \
4411  intermede *= x[var]; \
4412  } \
4413  continue; \
4414  }
4415  // case kpol:
4416  case kxpol: R__POLY(0);
4417  case kypol: R__POLY(1);
4418  case kzpol: R__POLY(2);
4419 
4420  case kDefinedVariable : {
4421  if (!precalculated) {
4422  precalculated = 1;
4423  for(j=0;j<fNval;j++) param_calc[j]=DefinedValue(j);
4424  }
4425  pos++; tab[pos-1] = param_calc[(oper & kTFOperMask)];
4426  continue;
4427  }
4428 
4429  case kDefinedString : {
4430  int param = (oper & kTFOperMask);
4431  if (!precalculated_str) {
4432  precalculated_str=1;
4433  for (j=0;j<fNstring;j++) string_calc[j]=DefinedString(j);
4434  }
4435  strpos++; stringStack[strpos-1] = string_calc[param];
4436  pos++; tab[pos-1] = 0;
4437  continue;
4438  }
4439 
4440  case kFunctionCall: {
4441  // an external function call
4442 
4443  int param = (oper & kTFOperMask);
4444  int fno = param / 1000;
4445  int nargs = param % 1000;
4446 
4447  // Retrieve the function
4448  TMethodCall *method = (TMethodCall*)fFunctions.At(fno);
4449 
4450  // Set the arguments
4451  method->ResetParam();
4452  if (nargs) {
4453  UInt_t argloc = pos-nargs;
4454  for(j=0;j<nargs;j++,argloc++,pos--) {
4455  method->SetParam(tab[argloc]);
4456  }
4457  }
4458  pos++;
4459  Double_t ret;
4460  method->Execute(ret);
4461  tab[pos-1] = ret; // check for the correct conversion!
4462 
4463  continue;
4464  };
4465  }
4466  if (!TestBit(kOptimizationError)) {
4467  SetBit(kOptimizationError);
4468  Warning("EvalParFast","Found an unsupported optmized opcode (%d)",oper >> kTFOperShift);
4469  }
4470  }
4471  Double_t result0 = tab[0];
4472  return result0;
4473 
4474 }
4475 
4476 ////////////////////////////////////////////////////////////////////////////////
4477 /// Pre compile function
4478 
4479 Int_t TFormula::PreCompile()
4480 {
4481  TString str = fTitle;
4482  if (str.Length()<3) return 1;
4483  if (str[str.Length()-1]!='+'&&str[str.Length()-2]!='+') return 1;
4484  str[str.Length()-2]=0;
4485  TString funName("preformula_");
4486  funName += fName;
4487  if (ROOT::v5::TFormulaPrimitive::FindFormula(funName)) return 0;
4488  TString fileName;
4489  fileName.Form("/tmp/%s.C",funName.Data());
4490 
4491  FILE *hf;
4492  hf = fopen(fileName.Data(),"w");
4493  if (hf == 0) {
4494  Error("PreCompile","Unable to open the file %s for writing.",fileName.Data());
4495  return 1;
4496  }
4497  fprintf(hf, "/////////////////////////////////////////////////////////////////////////\n");
4498  fprintf(hf, "// This code has been automatically generated \n");
4499  //
4500  fprintf(hf, "Double_t %s(Double_t *x, Double_t *p){",funName.Data());
4501  fprintf(hf, "return (%s);\n}",str.Data());
4502 
4503  // fprintf("ROOT::v5::TFormulaPrimitive::AddFormula(new ROOT::v5::TFormulaPrimitive(\"%s::%s\",\"%s::%s\",(ROOT::v5::TFormulaPrimitive::GenFunc0)%s::%s));\n",
4504  // clname,method->GetName(),clname,method->GetName(),clname,method->GetName());
4505  fclose(hf);
4506 
4507  return 0;
4508 
4509 
4510 }
4511 
4512 ////////////////////////////////////////////////////////////////////////////////
4513 /// static function to set the maximum value of 3 parameters
4514 ///
4515 /// - maxop : maximum number of operations
4516 /// - maxpar : maximum number of parameters
4517 /// - maxconst : maximum number of constants
4518 ///
4519 /// None of these parameters cannot be less than 10 (default is 1000)
4520 /// call this function to increase one or all maxima when processing
4521 /// very complex formula, eg TFormula::SetMaxima(100000,1000,1000000);
4522 /// If you process many functions with a small number of operations/parameters
4523 /// you may gain some memory and performance by decreasing these values.
4524 
4525 void TFormula::SetMaxima(Int_t maxop, Int_t maxpar, Int_t maxconst)
4526 {
4527  gMAXOP = TMath::Max(10,maxop);
4528  gMAXPAR = TMath::Max(10,maxpar);
4529  gMAXCONST = TMath::Max(10,maxconst);
4530 }
4531 
4532 ////////////////////////////////////////////////////////////////////////////////
4533 /// static function to get the maximum value of 3 parameters
4534 /// -maxop : maximum number of operations
4535 /// -maxpar : maximum number of parameters
4536 /// -maxconst : maximum number of constants
4537 
4538 void TFormula::GetMaxima(Int_t& maxop, Int_t& maxpar, Int_t& maxconst)
4539 {
4540  maxop = gMAXOP;
4541  maxpar = gMAXPAR;
4542  maxconst = gMAXCONST;
4543 }
4544 
4545  } // end namespace v5
4546 
4547 } // end namespace ROOT