source: trunk/MagicSoft/Mars/mdata/MDataChain.cc@ 5899

Last change on this file since 5899 was 5832, checked in by tbretz, 20 years ago
*** empty log message ***
File size: 22.6 KB
Line 
1/* ======================================================================== *\
2!
3! *
4! * This file is part of MARS, the MAGIC Analysis and Reconstruction
5! * Software. It is distributed to you in the hope that it can be a useful
6! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
7! * It is distributed WITHOUT ANY WARRANTY.
8! *
9! * Permission to use, copy, modify and distribute this software and its
10! * documentation for any purpose is hereby granted without fee,
11! * provided that the above copyright notice appear in all copies and
12! * that both that copyright notice and this permission notice appear
13! * in supporting documentation. It is provided "as is" without express
14! * or implied warranty.
15! *
16!
17!
18! Author(s): Thomas Bretz 04/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
19!
20! Copyright: MAGIC Software Development, 2000-2005
21!
22!
23\* ======================================================================== */
24
25/////////////////////////////////////////////////////////////////////////////
26//
27// MDataChain
28// ==========
29//
30// With this chain you can concatenate simple mathematical operations on
31// members of mars containers.
32//
33//
34// Rules
35// -----
36//
37// In the constructor you can give rule, like
38// "HillasSource.fDist / MHillas.fLength"
39// Where MHillas/HillasSource is the name of the parameter container in
40// the parameter list and fDist/fLength is the name of the data members
41// in the containers. The result will be fDist divided by fLength.
42//
43// In case you want to access a data-member which is a data member object
44// you can acces it with (Remark: it must derive from MParContainer):
45// "MCameraLV.fPowerSupplyA.fVoltagePos5V"
46//
47// You can also use parantheses:
48// "HillasDource.fDist / (MHillas.fLength + MHillas.fWidth)"
49//
50//
51// Operators
52// ---------
53//
54// The allowed operations are: +, -, *, /, %, ^
55//
56// While a%b returns the floating point reminder of a/b.
57// While a^b returns a to the power of b
58//
59// Warning: There is no priority rule build in. So better use parantheses
60// to get correct results. The rule is parsed/evaluated from the left
61// to the right, which means:
62//
63// "MHillas.fWidth + MHillas.fLength / HillasSource.fDist"
64//
65// is parses as
66//
67// "(MHillas.fWidth + MHillas.fLength) / HillasSource.fDist"
68//
69// You can also use mathmatical operators, eg:
70// "5*log10(MMcEvt.fEnergy*MHillas.fSize)"
71//
72// The allowed operators are:
73// exp(x) e^x
74// log(x) natural logarithm of x
75// pow10(x) 10^x
76// log2(x) logarithm of x to base two
77// log10(x) logarithm of x to base ten
78// cos(x) cosine of x
79// sin(x) sine of x
80// tan(x) tangent of x
81// cosh(x) hyperbolic cosine of x
82// sinh(x) hyperbolic sine of x
83// tanh(x) hyperbolic tangent of x
84// acosh(x) arc hyperbolic cosine of x
85// asinh(x) arc hyperbolic sine of x
86// atanh(x) arc hyperbolic tangent of x
87// acos(x) arc cosine (inverse cosine) of x
88// asin(x) arc sine (inverse sine) of x
89// atan(x) arc tangent (inverse tangent) of x
90// sqrt(x) square root of x
91// sqr(x) square of x
92// abs(x) absolute value of x, |x|
93// floor(x) round down to the nearest integer (floor(9.9)=9)
94// ceil(x) round up to the nearest integer (floor(9.1)=10)
95// round(x) round to the nearest integer
96// r2d(x) transform radians to degrees
97// d2r(x) transform degrees to radians
98// rand(x) returns a uniform deviate on the interval ( 0, x ].
99// (gRandom->Uniform(x) is returned)
100// randp(x) returns gRandom->Poisson(x)
101// rande(x) returns gRandom->Exp(x)
102// randi(x) returns gRandom->Integer(x)
103// randg(x) returns gRandom->Gaus(0, x)
104// randl(x) returns gRandom->Landau(0, x)
105// isnan(x) return 1 if x is NaN (Not a Number) otherwise 0
106// finite(x) return 1 if the number is a valid double (not NaN, inf)
107//
108// NaN (Not a Number) means normally a number which is to small to be
109// stored in a floating point variable (eg. 0<x<1e-56 or similar) or
110// a number which function is not defined (like asin(1.5))
111//
112// inf is the symbol for an infinite number.
113//
114//
115// Constants
116// ---------
117//
118// Constants are implemented in ParseDataMember, namely:
119// kPi: TMath::Pi()
120// kRad2Deg: 180/kPi
121// kDeg2Rad: kPi/180
122//
123// You can also defined constants which are defined in TMath by:
124// kLn10 for static Double_t TMath::Ln10();
125// kLogE for static Double_t TMath::LogE();
126// kRadToDeg for static Double_t TMath::RadToDeg();
127// kDegToRad for static Double_t TMath::DegToRad();
128// ...
129//
130// Remark:
131// In older root versions only Pi() and E() are implemented
132// in TMath.
133//
134//
135// Variable Parameters
136// ------------------------
137// If you want to use variables, eg for fits you can use [0], [1], ...
138// These values are initialized with 0 and set by calling
139// SetVariables().
140//
141//
142// Multi-argument functions
143// ------------------------
144// You can use multi-argument functions, too. The access is implemented
145// via TFormula, which slows down thing a lot. If you can avoid usage
146// of such expression you accelerate the access a lot. Example:
147// "TMath::Hypot(MHillas.fMeanX, MHillas.MeanY)"
148// "pow(MHillas.fMeanX*MHillas.MeanY, -1.2)"
149// It should be possible to use all functions which are accessible
150// via the root-dictionary.
151//
152//
153// REMARK:
154// - All the random functions are returning 0 if gRandom==0
155// - You may get better results if you are using a TRandom3
156//
157// To Do:
158// - The possibility to use other objects inheriting from MData
159// is missing.
160// - By automatic pre-adding parantheses to the rule it would be possible
161// to implement priorities for operators.
162//
163/////////////////////////////////////////////////////////////////////////////
164
165#include "MDataChain.h"
166
167#include <ctype.h> // isalnum, ...
168#include <stdlib.h> // strtod, ...
169
170#include <TRandom.h>
171#include <TMethodCall.h>
172
173#include "MLog.h"
174#include "MLogManip.h"
175
176#include "MDataList.h"
177#include "MDataValue.h"
178#include "MDataMember.h"
179#include "MDataFormula.h"
180#include "MDataElement.h"
181
182ClassImp(MDataChain);
183
184using namespace std;
185
186// --------------------------------------------------------------------------
187//
188// Constructor which takes a rule and a surrounding operator as argument
189//
190MDataChain::MDataChain(const char *rule, OperatorType_t op)
191 : fOperatorType(op)
192{
193 fName = "MDataChain";
194 fTitle = rule;
195
196 fMember=ParseString(rule, 1);
197}
198
199// --------------------------------------------------------------------------
200//
201// Constructor taking a rule as an argument. For more details see
202// class description
203//
204MDataChain::MDataChain(const char *rule, const char *name, const char *title)
205 : fMember(NULL), fOperatorType(kENoop)
206{
207 fName = name ? name : "MDataChain";
208 fTitle = title ? title : rule;
209
210 if (TString(rule).IsNull())
211 return;
212
213 *fLog << inf << "Parsing rule... " << flush;
214 if (!(fMember=ParseString(rule, 1)))
215 {
216 *fLog << err << dbginf << "Parsing '" << rule << "' failed." << endl;
217 return;
218 }
219 *fLog << inf << "found: " << GetRule() << endl;
220}
221
222// --------------------------------------------------------------------------
223//
224// PreProcesses all members in the list
225//
226Bool_t MDataChain::PreProcess(const MParList *pList)
227{
228 return fMember ? fMember->PreProcess(pList) : kFALSE;
229}
230
231// --------------------------------------------------------------------------
232//
233// Checks whether at least one member has the ready-to-save flag.
234//
235Bool_t MDataChain::IsReadyToSave() const
236{
237 *fLog << all << "fM=" << fMember << "/" << (int)fMember->IsReadyToSave() << " " << endl;
238 return fMember ? fMember->IsReadyToSave() : kFALSE;
239}
240
241// --------------------------------------------------------------------------
242//
243// Destructor. Delete filters.
244//
245MDataChain::~MDataChain()
246{
247 if (fMember)
248 delete fMember;
249}
250
251// --------------------------------------------------------------------------
252//
253// Returns the number of alphanumeric characters (including '.' and ';')
254// in the given string
255//
256Int_t MDataChain::IsAlNum(TString txt)
257{
258 int l = txt.Length();
259 for (int i=0; i<l; i++)
260 {
261 if (!isalnum(txt[i]) && txt[i]!='.' && txt[i]!=':' && txt[i]!=';' &&
262 /*txt[i]!='['&&txt[i]!=']'&&*/
263 ((txt[i]!='-' && txt[i]!='+') || i!=0))
264 return i;
265 }
266
267 return l;
268}
269
270Int_t MDataChain::GetBracket(TString txt, char open, char close)
271{
272 Int_t first=1;
273 for (int cnt=0; first<txt.Length(); first++)
274 {
275 if (txt[first]==open)
276 cnt++;
277 if (txt[first]==close)
278 cnt--;
279 if (cnt==-1)
280 break;
281 }
282 return first;
283}
284
285// --------------------------------------------------------------------------
286//
287// Compare a given text with the known operators. If no operator match
288// kENoop is retunred, otherwise the corresponding OperatorType_t
289//
290MDataChain::OperatorType_t MDataChain::ParseOperator(TString txt) const
291{
292 txt.ToLower();
293
294 if (txt=="abs") return kEAbs;
295 if (txt=="fabs") return kEAbs;
296 if (txt=="log") return kELog;
297 if (txt=="log2") return kELog2;
298 if (txt=="log10") return kELog10;
299 if (txt=="sin") return kESin;
300 if (txt=="cos") return kECos;
301 if (txt=="tan") return kETan;
302 if (txt=="sinh") return kESinH;
303 if (txt=="cosh") return kECosH;
304 if (txt=="tanh") return kETanH;
305 if (txt=="asin") return kEASin;
306 if (txt=="acos") return kEACos;
307 if (txt=="atan") return kEATan;
308 if (txt=="asinh") return kEASinH;
309 if (txt=="acosh") return kEACosH;
310 if (txt=="atanh") return kEATanH;
311 if (txt=="sqrt") return kESqrt;
312 if (txt=="sqr") return kESqr;
313 if (txt=="exp") return kEExp;
314 if (txt=="pow10") return kEPow10;
315 if (txt=="sgn") return kESgn;
316 if (txt=="floor") return kEFloor;
317 if (txt=="r2d") return kERad2Deg;
318 if (txt=="d2r") return kEDeg2Rad;
319 if (txt=="rand") return kERandom;
320 if (txt=="randp") return kERandomP;
321 if (txt=="rande") return kERandomE;
322 if (txt=="randi") return kERandomI;
323 if (txt=="randg") return kERandomG;
324 if (txt=="randl") return kERandomL;
325 if (txt=="isnan") return kEIsNaN;
326 if (txt=="finite") return kEFinite;
327 if (txt[0]=='-') return kENegative;
328 if (txt[0]=='+') return kEPositive;
329
330 return kENoop;
331}
332
333// --------------------------------------------------------------------------
334//
335// Here the names of data members are interpreted. This can be used to
336// check for constants.
337//
338MData *MDataChain::ParseDataMember(TString txt)
339{
340 //txt.ToLower();
341
342 if (txt=="kRad2Deg") return new MDataValue(kRad2Deg);
343 if (txt=="kDeg2Rad") return new MDataValue(1./kRad2Deg);
344
345 if (!txt.BeginsWith("k"))
346 return new MDataMember(txt.Data());
347
348 const TString name = txt(1, txt.Length());
349#if ROOT_VERSION_CODE < ROOT_VERSION(4,02,00)
350 TMethodCall call(TMath::Class(), name, "");
351#else
352 static TClass *const tmath = TClass::GetClass("TMath");
353 TMethodCall call(tmath, name, "");
354#endif
355
356 switch (call.ReturnType())
357 {
358 case TMethodCall::kLong:
359 Long_t l;
360 call.Execute(l);
361 return new MDataValue(l);
362 case TMethodCall::kDouble:
363 Double_t d;
364 call.Execute(d);
365 return new MDataValue(d);
366 default:
367 break;
368 }
369
370 return new MDataMember(txt.Data());
371}
372
373void MDataChain::SimplifyString(TString &txt) const
374{
375 while (txt.First("--")>=0 || txt.First("++")>=0 ||
376 txt.First("+-")>=0 || txt.First("-+")>=0)
377 {
378 txt.ReplaceAll("--", "+");
379 txt.ReplaceAll("++", "+");
380 txt.ReplaceAll("-+", "-");
381 txt.ReplaceAll("+-", "-");
382 }
383}
384
385// --------------------------------------------------------------------------
386//
387// Core of the data chain. Here the chain is constructed out of the rule.
388//
389MData *MDataChain::ParseString(TString txt, Int_t level)
390{
391 if (level==0)
392 SimplifyString(txt);
393
394 MData *member0=NULL;
395
396 char type=0;
397 int nlist = 0;
398
399 while (!txt.IsNull())
400 {
401 MData *newmember = NULL;
402
403 txt = txt.Strip(TString::kBoth);
404
405 switch (txt[0])
406 {
407 case '(':
408 {
409 //
410 // Search for the corresponding parantheses
411 //
412 const Int_t first=GetBracket(txt, '(', ')');
413 if (first==txt.Length())
414 {
415 *fLog << err << dbginf << "Syntax Error: ')' missing." << endl;
416 if (member0)
417 delete member0;
418 return NULL;
419 }
420
421 //
422 // Make a copy of the 'interieur' and delete the substring
423 // including the brackets
424 //
425 TString sub = txt(1, first-1);
426 txt.Remove(0, first+1);
427
428 //
429 // Parse the substring
430 //
431 newmember = ParseString(sub, level+1);
432 if (!newmember)
433 {
434 *fLog << err << dbginf << "Parsing '" << sub << "' failed." << endl;
435 if (member0)
436 delete member0;
437 return NULL;
438 }
439 }
440 break;
441
442 case ')':
443 *fLog << err << dbginf << "Syntax Error: Too many ')'" << endl;
444 if (member0)
445 delete member0;
446 return NULL;
447
448 case '+':
449 case '-':
450 case '*':
451 case '/':
452 case '%':
453 case '^':
454 if (member0)
455 {
456 //
457 // Check for the type of the symbol
458 //
459 char is = txt[0];
460 txt.Remove(0, 1);
461
462 //
463 // If no filter is available or the available filter
464 // is of a different symbols we have to create a new
465 // data list with the new symbol
466 //
467 if (/*!member0->InheritsFrom(MDataMember::Class()) ||*/ type!=is)
468 {
469 MDataList *list = new MDataList(is);
470 list->SetName(Form("List_%c_%d", is, 10*level+nlist++));
471
472 list->SetOwner();
473 list->AddToList(member0);
474
475 member0 = list;
476
477 type = is;
478 }
479 continue;
480 }
481
482 if (txt[0]!='-' && txt[0]!='+')
483 {
484 *fLog << err << dbginf << "Syntax Error: First argument of '";
485 *fLog << txt[0] << "' opartor missing." << endl;
486 if (member0)
487 delete member0;
488 return NULL;
489 }
490
491 // FALLTHROU
492
493 case '0':
494 case '1':
495 case '2':
496 case '3':
497 case '4':
498 case '5':
499 case '6':
500 case '7':
501 case '8':
502 case '9':
503 case '[':
504 case ']':
505 if ((txt[0]!='-' && txt[0]!='+') || isdigit(txt[1]) || txt[1]=='.')
506 {
507 if (!txt.IsNull() && txt[0]=='[')
508 {
509 Int_t first = GetBracket(txt, '[', ']');
510 TString op = txt(1, first-1);
511 txt.Remove(0, first+1);
512
513 newmember = new MDataValue(0, atoi(op));
514 break;
515 }
516
517 char *end;
518 Double_t num = strtod(txt.Data(), &end);
519 if (!end || txt.Data()==end)
520 {
521 *fLog << err << dbginf << "Error trying to convert '" << txt << "' to value." << endl;
522 if (member0)
523 delete member0;
524 return NULL;
525 }
526
527 txt.Remove(0, end-txt.Data());
528
529 newmember = new MDataValue(num);
530 break;
531 }
532
533 // FALLTHROUH
534
535 default:
536 int i = IsAlNum(txt);
537
538 if (i==0)
539 {
540 *fLog << err << dbginf << "Syntax Error: Name of data member missing in '" << txt << "'" << endl;
541 if (member0)
542 delete member0;
543 return NULL;
544 }
545
546 TString text = txt(0, i);
547
548 txt.Remove(0, i);
549 txt = txt.Strip(TString::kBoth);
550
551 if (!txt.IsNull() && txt[0]=='[')
552 {
553 Int_t first = GetBracket(txt, '[', ']');
554 TString op = txt(1, first-1);
555 txt.Remove(0, first+1);
556
557 newmember = new MDataElement(text, atoi(op));
558 break;
559 }
560 if ((txt.IsNull() || txt[0]!='(') && text[0]!='-' && text[0]!='+')
561 {
562 newmember = ParseDataMember(text.Data());
563 break;
564 }
565
566 OperatorType_t op = ParseOperator(text);
567
568 Int_t first = GetBracket(txt, '(', ')');
569 TString sub = op==kENegative || op==kEPositive ? text.Remove(0,1) + txt : txt(1, first-1);
570 txt.Remove(0, first+1);
571
572 if (op==kENoop)
573 {
574 newmember = new MDataFormula(Form("%s(%s)", (const char*)text, (const char*)sub));
575 if (newmember->IsValid())
576 break;
577
578 *fLog << err << dbginf << "Syntax Error: Operator '" << text << "' unknown." << endl;
579 if (member0)
580 delete member0;
581 return NULL;
582 }
583
584 newmember = new MDataChain(sub, op);
585 if (!newmember->IsValid())
586 {
587 *fLog << err << dbginf << "Syntax Error: Error parsing contents '" << sub << "' of operator " << text << endl;
588 delete newmember;
589 if (member0)
590 delete member0;
591 return NULL;
592 }
593 }
594
595 if (!member0)
596 {
597 member0 = newmember;
598 continue;
599 }
600
601 if (!member0->InheritsFrom(MDataList::Class()))
602 continue;
603
604 ((MDataList*)member0)->AddToList(newmember);
605 }
606
607 return member0;
608}
609
610// --------------------------------------------------------------------------
611//
612// Returns the value described by the rule
613//
614Double_t MDataChain::GetValue() const
615{
616 if (!fMember)
617 {
618 //*fLog << warn << "MDataChain not valid." << endl;
619 return 0;
620 }
621
622 const Double_t val = fMember->GetValue();
623
624 switch (fOperatorType)
625 {
626 case kEAbs: return TMath::Abs(val);
627 case kELog: return TMath::Log(val);
628 case kELog2: return TMath::Log2(val);
629 case kELog10: return TMath::Log10(val);
630 case kESin: return TMath::Sin(val);
631 case kECos: return TMath::Cos(val);
632 case kETan: return TMath::Tan(val);
633 case kESinH: return TMath::SinH(val);
634 case kECosH: return TMath::CosH(val);
635 case kETanH: return TMath::TanH(val);
636 case kEASin: return TMath::ASin(val);
637 case kEACos: return TMath::ACos(val);
638 case kEATan: return TMath::ATan(val);
639 case kEASinH: return TMath::ASinH(val);
640 case kEACosH: return TMath::ACosH(val);
641 case kEATanH: return TMath::ATanH(val);
642 case kESqrt: return TMath::Sqrt(val);
643 case kESqr: return val*val;
644 case kEExp: return TMath::Exp(val);
645 case kEPow10: return TMath::Power(10, val);
646 case kESgn: return val<0 ? -1 : 1;
647 case kENegative: return -val;
648 case kEPositive: return val;
649 case kEFloor: return TMath::Floor(val);
650 case kECeil: return TMath::Ceil(val);
651 case kERound: return TMath::Nint(val);
652 case kERad2Deg: return val*180/TMath::Pi();
653 case kEDeg2Rad: return val*TMath::Pi()/180;
654 case kERandom: return gRandom ? gRandom->Uniform(val) : 0;
655 case kERandomP: return gRandom ? gRandom->Poisson(val) : 0;
656 case kERandomE: return gRandom ? gRandom->Exp(val) : 0;
657 case kERandomI: return gRandom ? gRandom->Integer((int)val) : 0;
658 case kERandomG: return gRandom ? gRandom->Gaus(0, val) : 0;
659 case kERandomL: return gRandom ? gRandom->Landau(0, val) : 0;
660 case kEIsNaN: return TMath::IsNaN(val);
661 case kEFinite: return TMath::Finite(val);
662 case kENoop: return val;
663 }
664
665 *fLog << warn << "No Case for " << fOperatorType << " available." << endl;
666
667 return 0;
668}
669
670// --------------------------------------------------------------------------
671//
672// Builds a rule from all the chain members. This is a rule which could
673// be used to rebuild the chain.
674//
675TString MDataChain::GetRule() const
676{
677 if (!fMember)
678 return "<n/a>";
679
680 TString str;
681
682 Bool_t bracket = fOperatorType!=kENoop && !fMember->InheritsFrom(MDataList::Class());
683
684 switch (fOperatorType)
685 {
686 case kEAbs: str += "abs" ; break;
687 case kELog: str += "log" ; break;
688 case kELog2: str += "log2" ; break;
689 case kELog10: str += "log10" ; break;
690 case kESin: str += "sin" ; break;
691 case kECos: str += "cos" ; break;
692 case kETan: str += "tan" ; break;
693 case kESinH: str += "sinh" ; break;
694 case kECosH: str += "cosh" ; break;
695 case kETanH: str += "tanh" ; break;
696 case kEASin: str += "asin" ; break;
697 case kEACos: str += "acos" ; break;
698 case kEATan: str += "atan" ; break;
699 case kEASinH: str += "asinh" ; break;
700 case kEACosH: str += "acosh" ; break;
701 case kEATanH: str += "atanh" ; break;
702 case kESqrt: str += "sqrt" ; break;
703 case kESqr: str += "sqr" ; break;
704 case kEExp: str += "exp" ; break;
705 case kEPow10: str += "pow10" ; break;
706 case kESgn: str += "sgn" ; break;
707 case kENegative: str += "-" ; break;
708 case kEPositive: str += "+" ; break;
709 case kEFloor: str += "floor" ; break;
710 case kECeil: str += "ceil" ; break;
711 case kERound: str += "round" ; break;
712 case kERad2Deg: str += "r2d" ; break;
713 case kEDeg2Rad: str += "d2r" ; break;
714 case kERandom: str += "rand" ; break;
715 case kERandomP: str += "randp" ; break;
716 case kERandomE: str += "rande" ; break;
717 case kERandomI: str += "randi" ; break;
718 case kERandomG: str += "randg" ; break;
719 case kERandomL: str += "randl" ; break;
720 case kEIsNaN: str += "isnan" ; break;
721 case kEFinite: str += "finite"; break;
722 case kENoop:
723 break;
724 }
725
726 if (bracket)
727 str += "(";
728
729 str += fMember->GetRule();
730
731 if (bracket)
732 str += ")";
733
734 return str;
735}
736
737// --------------------------------------------------------------------------
738//
739// Return a comma seperated list of all data members used in the chain.
740// This is mainly used in MTask::AddToBranchList
741//
742TString MDataChain::GetDataMember() const
743{
744 return fMember->GetDataMember();
745}
746
747void MDataChain::SetVariables(const TArrayD &arr)
748{
749 return fMember->SetVariables(arr);
750}
Note: See TracBrowser for help on using the repository browser.