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

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