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

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